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FOREWORD 


It  has  been  known  for  20  years  that  LES  predicts  incorrectly  mean  velocity  gradients  near  the  ground.  This  serious  error  causes 
inaccuracy  in  other  important  predictions  including,  in  the  presence  of  convection,  the  entire  structure  of  the  ABL,  heat,  humidity 
and  contaminant  transport,  and  acoustic  propagation.  Many  studies  have  tried  to  eliminate  the  overshoot  through  the  SFS 
model,  but  none  have  been  entirely  successful  because  the  mechanisms  are  not  understood.  We  have  finally  succeeded  in 
explaining  the  mechanisms  underlying  the  overshoot  and  we  have  defined  the  requirements  for  high-accuracy  LES  of  the  ABL 
surface  layer.  These  requirements  include  a  mix  of  physical,  numerical  and  modeling  issues.  Underlying  the  explanation  is  a 
crucial  distinction  between  “numerical  friction”  and  real  friction  and  their  similar  roles.  We  show  that  numerical  friction  includes  a 
mix  of  modeling,  grid,  and  numerical  elements  that  underlie  the  overshoot  as  well  as  poor  predictions  of  law-of-the-wall.  The 
interaction  among  the  elements  underlying  numerical  friction  underlies  the  overshoot  and  high-accuracy  LES  of  the  surface 
layer.  Underlying  these  primary  issues  are  related  issues  that  must  be  addressed  in  high-accuracy  LES,  including  numerical 
instability,  lower  boundary  condition,  and  details  of  the  SFS  closure  prediction  and  grid. 

STATEMENT  OF  THE  PROBLEM  STUDIED 

In  1992  Mason  and  Thomson  made  the  critical  observation  that  their  large-eddy  simulations  (LES)  of  the  atmospheric  boundary 
layer  (ABL)  did  not  predict  correctly  the  mean  velocity  profile  near  the  ground.  They  observed  that  gradients  of  mean  velocity 
are  much  higher  than  predicted  by  law-of-the-wall  scaling  and  shown  by  data  in  the  part  of  the  computational  domain  adjacent 
to  the  ground.  They  assumed  that  the  over-prediction  was  a  consequence  of  the  failure  of  the  Smagorinsky  model,  and  showed 
that  the  over-prediction  could  be  altered  (but  not  removed)  by  adding  a  random  solenoidal  acceleration  to  the  Navier-Stokes 
equation,  in  effect  randomizing  the  subfilter-scale  (SFS)  stress  divergence.  In  the  intervening  years  it  has  been  discovered  that 
this  over-prediction  is  a  fundamental  flaw  in  LES  of  the  ABL  when  the  surface  layer  is  shear-dominated,  and  is  particularly 
associated  with  eddy  viscosity  representations  of  SFS  stress.  A  number  of  researchers  have  found  that  it  is  possible  to  reduce 
the  degree  of  error  either  by  adjusting  details  of  the  eddy  viscosity  representation  or  by  adjusting  the  eddy  viscosity  closure 
itself.  However,  none  have  succeeded  in  fully  eliminating  the  error,  primarily  because  the  mechanisms  underlying  the  flaw  are 
not  understood. 

The  issue  is  much  more  serious  in  predicting  geophysical  boundary  layers  than  simply  an  inaccurate  prediction  of  mean  velocity 
near  the  ground.  Over-prediction  of  mean  shear  affects  the  prediction  of  turbulent  kinetic  energy  production  and  directly  alters 
turbulence  structure  and  turbulence-driven  fluxes  in  the  part  of  the  planetary  boundary  layer  in  which  humans  reside. 
Consequently,  the  application  of  LES  to  problems  that  require  accurate  prediction  of  meteorological  events  near  the  ground  is 
questionable.  Examples  include  ground-to-ground  acoustic  propagation,  dispersion  of  contaminants  originating  at  the  ground, 
and  wind  speed  and  direction.  As  importantly,  it  has  been  shown  that  because  atmospheric  thermals  originate  at  the  ground  in 
the  part  of  the  boundary  layer  where  mean  shear  is  over-predicted,  the  LES  predictions  of  thermal  structure  are  also  incorrect 
under  moderately  convective  conditions,  not  only  near  the  surface  but  throughout  the  ABL,  and  that  these  errors  are  a  direct 
consequence  of  the  over-prediction  of  mean  shear  at  the  surface.  Consequently,  LES  predictions  of  transport  from  the  lower  to 
upper  ABL — humidity,  heat,  contaminants  and  other  scalars — are  fundamentally  inaccurate,  and  could  adversely  affect  other 
predictions  such  as  cloud  formation  and  long-range  scalar  dispersion. 

Given  the  importance  of  near-surface  predictions  in  geophysical  flows,  our  program  of  research  focused  on  resolving  the  now 
23-year-old  problem  of  inaccurate  LES  prediction  of  mean  gradients  in  the  surface  layer. 

SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 

After  a  long  period  of  methodical  detective  work  and  a  combination  of  mathematical  and  physical  analysis  combined  with  an 
extensive  large-eddy  simulation  (LES)  campaign  of  atmospheric  boundary  layer  (ABL)  simulations,  we  succeeded  in  solving  the 
puzzle  of  why  it  is  that  large-eddy  simulation  of  the  atmospheric  boundary  layer  over-predicts  mean  shear  adjacent  to  the 
ground,  and  we  have  learned  what  elements  are  needed  to  address  the  problem  and  improve  accuracy  of  surface  layer 
predictions.  As  might  be  expected,  the  solution  is  multi-faceted,  and  we  developed  a  series  of  requirements  for  high-accuracy 
LES  of  the  ABL  that  involve  a  mix  of  numerical,  gridding,  and  modeling  issues.  From  our  analysis  we  order  the  solution  process 
into  a  hierarchy  of  inter-related  issues  that  must  be  addressed  for  high-accuracy  LES  of  the  surface  layer. 

A  highly  comprehensive  and  though  explanation  of  the  background  to  the  fundamental  problem  in  LES,  the  mathematical  theory 
we  developed,  the  "R-ReLES"  parameter  space  we  created,  and  the  analysis  of  over  100  ABL  LES  simulations  to  place  results 
within  this  parameter  space  and  validate  both  the  theory  and  the  fundamental  explanation  for  the  error,  is  given  in  an  invited 
paper  by  the  Physics  of  Fluids  entitled  "Designing  Large-Eddy  Simulation  of  the  Turbulent  Boundary  Layer  to  Capture  Law-of- 
the-Wall  Scaling"  by  Brasseur  &  Wei.  The  paper  appeared  in  2010  and  has  attracted  a  great  deal  of  attention  within  the  large 
community  that  applies  LES  to  the  simulation  of  high  Reynolds  number  wall  bounded  flows.  This  work  has  also  lead  to  a  large 
number  of  invited  presentations  at  universities,  national  labs  and  meetings. 


Deviations  from  the  law-of-the-wall  in  mean  velocity  and  velocity  gradient  in  the  inertial  surface  layer  arise  from  a  competition 
between  the  characteristic  velocity  and  length  scales  and  other  velocity  or  length  scales  that  enter  either  from  the  true  fluid 
physics  or  during  the  conversion  to  the  discretized  dynamical  system  that  is  ultimately  advanced  on  the  computer.  Turbulence 
motions  in  the  surface  layer  at  the  characteristic  length  scale  z,  for  example,  may  compete  with  motions  at  the  boundary  layer 
scale  (Khanna  &  Brasseur  1997)  or  with  the  influences  of  surface  friction  or  roughness  that  are  characterized  by  the  viscous 
and  roughness  length  scales.  If  sufficiently  strong,  these  confounding  scales  will  alter  the  scaling  of  the  mean  velocity  gradient 
in  the  inertial  surface  layer.  Law-of-the-wall  assumes  that  this  is  not  the  case  and  experiments  in  the  laboratory  and  atmosphere 
have  generally  supported  this  assumption  when  there  are  no  confounding  scales  and  the  ratio  of  the  outer  to  viscous  length  and 
roughness  scales  are  sufficiently  large.  However,  when  the  influence  of  confounding  characteristic  scales  is  sufficiently  strong 
to  compete  with  the  characteristic  velocity  and  length  scales  in  the  inertial  surface  layer,  deviations  from  LOTW  result. 

What  we  have  shown  is  that,  in  the  conversion  from  the  exact  continuous  dynamical  system  to  the  LES  discretized  dynamical 
system  that  is  actually  advanced  on  the  computer,  additional  characteristic  scales  are  introduced  into  the  simulated  dynamics 
that  can  interfere  with  LOTW  scaling  in  the  surface  layer  when  sufficiently  strong.  Several  elements  in  the  simulation  might 
introduce  spurious  characteristic  scales:  (i)  models  of  existing  terms  and  new  terms  that  are  introduced  into  the  governing 
equation,  (ii)  models  for  unknown  boundary  conditions,  (iii)  the  type  and  order  of  the  discretization  of  derivatives  together  with 
grid  geometry,  and  (vi)  algorithmic  additions,  for  example,  to  maintain  numerical  stability.  The  current  study  has  focused 
primarily  on  the  influences  of  a  spurious  viscous  length  scale  that  is  introduced  within  the  computational  domain  by  the  model 
for  the  SFS  stress  tensor.  This  spurious  scale  is  a  reflection  of  the  manner  in  which  the  net  transfer  of  turbulence  energy  from 
resolved  to  subfilter  scales  is  modeled.  Whereas  the  interscale  interactions  that  underlie  the  transfer  of  energy  are  purely 
inertial  in  reality,  all  practical  closures  for  the  SFS  stress  tensor  model  the  net  transfer  of  energy  from  the  resolved  scales  by  a 
dissipative  mechanism  that  removes  energy  at  the  smallest  resolved  scales. 

However,  the  SFS  stress  is  often  not  the  only  term  that  is  modeled  in  the  conversion  from  the  continuous  to  the  discrete 
dynamical  system.  At  the  very  high  Reynolds  numbers  of  practical  interest,  and  at  which  law-of-the-wall  is  valid,  the  surface 
viscous  layer,  if  it  exists,  cannot  be  resolved  in  a  large-eddy  simulation.  The  vertical  derivatives  in  resolved  and  SFS  stress 
therefore  require  that  a  model  be  supplied  for  the  total  stress  at  the  surface.  Hybrid  schemes  couple  a  RANS  SFS  stress  model 
on  a  very  high  aspect  ratio  RANS  grid  in  an  extremely  thin  surface  viscous  layer  with  large-eddy  simulation  and  SFS  closure  on 
an  LES  grid  beginning  in  the  lower  inertial  surface  layer,  also  very  near  the  surface.  Thus,  the  LES  part  of  the  simulation  obtains 
the  lower  boundary  condition  on  stress  from  the  RANS  part  of  the  simulation,  which  is  far  from  exact.  Nonhybrid  schemes  must 
supply  a  model  for  the  total  horizontal  shear  stress,  generally  written  as  a  function  of  resolved  velocity  within  the  computational 
domain.  So,  in  addition  to  any  spurious  scales  introduced  within  the  computational  domain  by  the  SFS  stress  model,  there 
exists  the  possibility  that  the  model  for  surface  shear  stress  might  introduce  spurious  scales  at  the  boundary  of  the  LES 
computational  domain. 

We  have  found  experimentally  that  standard  surface  stress  models  do  introduce  spurious  effects  that  force  deviations  from 
LOTW  at  the  first  couple  grid  levels  adjacent  to  the  surface.  Fig.  10  shows  that  this  additional  confounding  contribution  from  the 
lower  stress  boundary  condition  is  obscured  when  the  frictional  contribution  from  the  interior  SFS  stress  is  sufficiently 
overwhelming  to  produce  the  overshoot.  When  the  LES  is  moved  into  the  HAZ  so  that  the  viscous  effects  causing  the  overshoot 
are  suppressed,  the  confounding  influences  of  the  surface  stress  model  become  apparent  and  spread  vertically  as  the  relative 
contribution  from  the  interior  SFS  stress  model  diminishes  with  increasing  R.  We  have  also  shown  that  adjustments  to  surface 
shear  stress  are  possible  which  significantly  reduce  the  spurious  influence  of  standard  surface  stress  models 
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FOREWORD 

It  has  been  known  for  20  years  that  large-eddy  simulation  (LES)  predicts  incorrectly  mean  velocity 
gradients  near  the  ground.  This  serious  error  causes  inaccuracy  in  other  important  predictions  including, 
in  the  presence  of  convection,  the  entire  structure  of  the  atmospheric  boundary  layer  (ABL),  heat, 
humidity  and  contaminant  transport,  and  acoustic  propagation.  Many  studies  have  tried  to  eliminate  the 
overshoot  through  the  SFS  model,  but  none  have  been  entirely  successful  because  the  mechanisms  are  not 
understood.  We  have  finally  succeeded  in  explaining  the  mechanisms  underlying  the  overshoot  and  we 
have  defined  the  requirements  for  high-accuracy  LES  of  the  ABL  surface  layer.  These  requirements 
include  a  mix  of  physical,  numerical  and  modeling  issues.  Underlying  the  explanation  is  a  crucial 
distinction  between  “numerical  friction”  and  real  friction  and  their  similar  roles.  We  show  that  numerical 
friction  includes  a  mix  of  modeling,  grid,  and  numerical  elements  that  underlie  the  overshoot  as  well  as 
poor  predictions  of  law-of-the-wall.  The  interaction  among  the  elements  underlying  numerical  friction 
underlies  the  overshoot  and  high-accuracy  LES  of  the  surface  layer.  Underlying  these  primary  issues  are 
related  issues  that  must  be  addressed  in  high-accuracy  LES,  including  numerical  instability,  lower 
boundary  condition,  and  details  of  the  SFS  closure  prediction  and  grid. 


STATEMENT  OF  THE  PROBLEM  STUDIED 


In  1992  Mason  and  Thomson  made  the  critical  observation  that  their  large-eddy  simulations  of  the 
atmospheric  boundary  layer  (ABL)  did  not  predict  correctly  the  mean  velocity  profile  near  the  ground. 
The  pointed  out  that  the  mean  velocity  gradient  scaled  according  to  inertial  Law-of-the-Wall  (LOTW) 
scaling, 


(/>(z)  = 


^dU_ 
u  *  dz 


(1) 


produces  a  well-defined  peak  within  the  surface  layer;  the  LES  does  not  predict  the  LOTW.  This 
“overshoot”  has  since  been  pointed  out  and  studied  by  many  researchers. 

Mason  and  Thomson  assumed  that  the  over  prediction  of  gradients  of  mean  velocity  relative  to  law- 
of-the-wall  scaling  and  shown  by  data  in  the  part  of  the  computational  domain  adjacent  to  the  ground  is  a 
consequence  of  the  failure  of  the  Smagorinsky  model,  and  showed  that  the  over -prediction  could  be 
altered  (but  not  removed)  by  adding  a  random  solenoidal  acceleration  to  the  Navier-Stokes  equation,  in 
effect  randomizing  the  subfilter-scale  (SFS)  stress  divergence.  In  the  intervening  years  it  has  been 
discovered  that  this  over-prediction  is  a  fundamental  flaw  in  LES  of  the  ABL  when  the  surface  layer  is 
shear-dominated,  and  is  particularly  associated  with  eddy  viscosity  representations  of  SFS  stress.  A 
number  of  researchers  have  found  that  it  is  possible  to  reduce  the  degree  of  error  either  by  adjusting 
details  of  the  eddy  viscosity  representation  or  by  adjusting  the  eddy  viscosity  closure  itself.  However, 
none  have  succeeded  in  fully  eliminating  the  error,  primarily  because  the  mechanisms  underlying  the  flaw 
are  not  understood. 

The  issue  is  much  more  serious  in  predicting  geophysical  boundary  layers  than  simply  an  inaccurate 
prediction  of  mean  velocity  near  the  ground.  Over-prediction  of  mean  shear  affects  the  prediction  of 
turbulent  kinetic  energy  production  and  directly  alters  turbulence  structure  and  turbulence-driven  fluxes  in 
the  part  of  the  planetary  boundary  layer  in  which  humans  reside.  Consequently,  the  application  of  LES  to 
problems  that  require  accurate  prediction  of  meteorological  events  near  the  ground  is  questionable. 
Examples  include  ground-to-ground  acoustic  propagation,  dispersion  of  contaminants  originating  at  the 
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ground,  and  wind  speed  and  direction.  As  importantly,  it  has  been  shown  that  because  atmospheric 
thermals  originate  at  the  ground  in  the  part  of  the  boundary  layer  where  mean  shear  is  over-predicted,  the 
LES  predictions  of  thermal  structure  are  also  incorrect  under  moderately  convective  conditions,  not  only 
near  the  surface  but  throughout  the  ABL,  and  that  these  errors  are  a  direct  consequence  of  the  over¬ 
prediction  of  mean  shear  at  the  surface.  Consequently,  LES  predictions  of  transport  from  the  lower  to 
upper  ABL — humidity,  heat,  contaminants  and  other  scalars — are  fundamentally  inaccurate,  and  could 
adversely  affect  other  predictions  such  as  cloud  formation  and  long-range  scalar  dispersion. 

Given  the  importance  of  near-surface  predictions  in  geophysical  flows,  our  program  of  research 
focused  on  resolving  the  now  23 -year-old  problem  of  inaccurate  LES  prediction  of  mean  gradients  in  the 
surface  layer. 

SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 

After  a  long  period  of  methodical  detective  work  and  a  combination  of  mathematical  and  physical 
analysis  combined  with  an  extensive  large-eddy  simulation  (LES)  campaign  of  atmospheric  boundary 
layer  (ABL)  simulations,  we  succeeded  in  solving  the  puzzle  of  why  it  is  that  large-eddy  simulation  of  the 
atmospheric  boundary  layer  over-predicts  mean  shear  adjacent  to  the  ground,  and  we  have  learned  what 
elements  are  needed  to  address  the  problem  and  improve  accuracy  of  surface  layer  predictions.  As  might 
be  expected,  the  solution  is  multi-faceted,  and  we  developed  a  series  of  requirements  for  high-accuracy 
LES  of  the  ABL  that  involve  a  mix  of  numerical,  gridding,  and  modeling  issues.  From  our  analysis  we 
order  the  solution  process  into  a  hierarchy  of  inter-related  issues  that  must  be  addressed  for  high-accuracy 
LES  of  the  surface  layer. 

A  highly  comprehensive  and  though  explanation  of  the  background  to  the  fundamental  problem  in 
LES,  the  mathematical  theory  we  developed,  the  ”R-ReLES”  parameter  space  we  created,  and  the 
analysis  of  over  100  ABL  LES  simulations  to  place  results  within  this  parameter  space  and  validate  both 
the  theory  and  the  fundamental  explanation  for  the  error,  is  given  in  an  invited  paper  by  the  Physics  of 
Fluids  entitled  "Designing  Large-Eddy  Simulation  of  the  Turbulent  Boundary  Layer  to  Capture  Law-of- 
the-Wall  Scaling”  by  Brasseur  &  Wei.  The  paper  appeared  in  2010  and  has  attracted  a  great  deal  of 
attention  within  the  large  community  that  applies  LES  to  the  simulation  of  high  Reynolds  number  wall 
bounded  flows.  This  work  has  also  lead  to  a  large  number  of  invited  presentations  at  universities,  national 
labs  and  meetings. 

Deviations  from  the  law-of-the-wall  in  mean  velocity  and  velocity  gradient  in  the  inertial  surface 
layer  arise  from  a  competition  between  the  characteristic  velocity  and  length  scales  and  other  velocity  or 
length  scales  that  enter  either  from  the  true  fluid  physics  or  during  the  conversion  to  the  discretized 
dynamical  system  that  is  ultimately  advanced  on  the  computer.  Turbulence  motions  in  the  surface  layer  at 
the  characteristic  length  scale  z,  for  example,  may  compete  with  motions  at  the  boundary  layer  scale 
(Khanna  &  Brasseur  1997)  or  with  the  influences  of  surface  friction  or  roughness  that  are  characterized 
by  the  viscous  and  roughness  length  scales.  If  sufficiently  strong,  these  confounding  scales  will  alter  the 
scaling  of  the  mean  velocity  gradient  in  the  inertial  surface  layer.  Law-of-the-wall  assumes  that  this  is  not 
the  case  and  experiments  in  the  laboratory  and  atmosphere  have  generally  supported  this  assumption 
when  there  are  no  confounding  scales  and  the  ratio  of  the  outer  to  viscous  length  and  roughness  scales  are 
sufficiently  large.  However,  when  the  influence  of  confounding  characteristic  scales  is  sufficiently  strong 
to  compete  with  the  characteristic  velocity  and  length  scales  in  the  inertial  surface  layer,  deviations  from 
LOTW  result. 

What  we  have  shown  is  that,  in  the  conversion  from  the  exact  continuous  dynamical  system  to  the 
LES  discretized  dynamical  system  that  is  actually  advanced  on  the  computer,  additional  characteristic 
scales  are  introduced  into  the  simulated  dynamics  that  can  interfere  with  LOTW  scaling  in  the  surface 
layer  when  sufficiently  strong.  Several  elements  in  the  simulation  might  introduce  spurious  characteristic 
scales:  (i)  models  of  existing  terms  and  new  terms  that  are  introduced  into  the  governing  equation,  (ii) 
models  for  unknown  boundary  conditions,  (iii)  the  type  and  order  of  the  discretization  of  derivatives 
together  with  grid  geometry,  and  (vi)  algorithmic  additions,  for  example,  to  maintain  numerical  stability. 
The  current  study  has  focused  primarily  on  the  influences  of  a  spurious  viscous  length  scale  that  is 
introduced  within  the  computational  domain  by  the  model  for  the  SFS  stress  tensor.  This  spurious  scale  is 
a  reflection  of  the  manner  in  which  the  net  transfer  of  turbulence  energy  from  resolved  to  subfilter  scales 
is  modeled.  Whereas  the  interscale  interactions  that  underlie  the  transfer  of  energy  are  purely  inertial  in 
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reality,  all  practical  closures  for  the  SFS  stress  tensor  model  the  net  transfer  of  energy  from  the  resolved 
scales  by  a  dissipative  mechanism  that  removes  energy  at  the  smallest  resolved  scales. 

However,  the  SFS  stress  is  often  not  the  only  term  that  is  modeled  in  the  conversion  from  the 
continuous  to  the  discrete  dynamical  system.  At  the  very  high  Reynolds  numbers  of  practical  interest,  and 
at  which  law-of-the-wall  is  valid,  the  surface  viscous  layer,  if  it  exists,  cannot  be  resolved  in  a  large-eddy 
simulation.  The  vertical  derivatives  in  resolved  and  SFS  stress  therefore  require  that  a  model  be  supplied 
for  the  total  stress  at  the  surface.  Hybrid  schemes  couple  a  RANS  SFS  stress  model  on  a  very  high  aspect 
ratio  RANS  grid  in  an  extremely  thin  surface  viscous  layer  with  large-eddy  simulation  and  SFS  closure 
on  an  LES  grid  beginning  in  the  lower  inertial  surface  layer,  also  very  near  the  surface.  Thus,  the  LES 
part  of  the  simulation  obtains  the  lower  boundary  condition  on  stress  from  the  RANS  part  of  the 
simulation,  which  is  far  from  exact.  Nonhybrid  schemes  must  supply  a  model  for  the  total  horizontal 
shear  stress,  generally  written  as  a  function  of  resolved  velocity  within  the  computational  domain.  So,  in 
addition  to  any  spurious  scales  introduced  within  the  computational  domain  by  the  SFS  stress  model, 
there  exists  the  possibility  that  the  model  for  surface  shear  stress  might  introduce  spurious  scales  at  the 
boundary  of  the  LES  computational  domain. 

We  have  found  experimentally  that  standard  surface  stress  models  do  introduce  spurious  effects  that 
force  deviations  from  LOTW  at  the  first  couple  grid  levels  adjacent  to  the  surface.  Fig.  10  shows  that  this 
additional  confounding  contribution  from  the  lower  stress  boundary  condition  is  obscured  when  the 
frictional  contribution  from  the  interior  SFS  stress  is  sufficiently  overwhelming  to  produce  the  overshoot. 
When  the  LES  is  moved  into  the  HAZ  so  that  the  viscous  effects  causing  the  overshoot  are  suppressed, 
the  confounding  influences  of  the  surface  stress  model  become  apparent  and  spread  vertically  as  the 
relative  contribution  from  the  interior  SFS  stress  model  diminishes  with  increasing  R.  We  have  also 
shown  that  adjustments  to  surface  shear  stress  are  possible  which  significantly  reduce  the  spurious 
influence  of  standard  surface  stress  models 


DETAILS  OF  THE  ANALYSES  DEVELOPED  AND  SOLUTIONS  ACHIEVED 
I.  Background  on  “the  Overshoot  problem” 

Before  developing  our  analysis  into  essential  mechanisms  underlying  the  prediction  of  LOTW  scaling 
of  mean  velocity  gradient  in  the  surface  layer  with  large-eddy  simulation,  we  set  the  stage  with  discussion 
into  the  overshoot  and  its  history.  Previous  studies,  many  of  which  have  produced  significant  reductions 
in  the  degree  of  overshoot,  have  provided  important  clues  that  lead  to  the  theory  developed  here,  and 
present  results  that  a  theory  should  explain. 

A.  Consequences  of  an  overshoot  in  mean  shear  rate 

Figure  1  shows  essential  aspects  of  the  overshoot  from  several  previous  large-eddy  simulations  of  the 
neutral  boundary  layer  that  have  focused  on  this  issue  since  Mason  &  Thomson1.  Because  the  overshoot 
is  associated  with  the  shear-dominated  region  of  the  boundary  layer,  it  is  particularly  apparent  in  the  fully 
shear-driven  neutral  boundary  layer  (Fig.  1).  Piomelli  and  Balaras16  point  out  that  in  DES  of  the  shear- 
driven  boundary  layer,  “unphysical,  nearly  one-dimensional,  wall  streaks  were  present  in  the  RANS 
region  . . .  and  shorter-scale  outer-layer  eddies  were  progressively  formed  as  one  moved  away  from  the 
wall.”  Similar  observations  have  been  made  in  the  neutral  atmospheric  boundary  layer17. 

In  the  moderately  convective  atmospheric  boundary  layer  an  overshoot  is  produced  in  the  shear- 
dominated  surface  layer  while  the  mixed  layer  is  buoyancy  dominated18.  In  the  presence  of  convection, 
vertically  driven  thermals  couple  the  surface  and  outer  boundary  layers.  Khanna  &  Brasseur17  showed  that 
the  elongated  structure  of  streamwise  turbulence  fluctuations  generated  near  the  surface  by  the  interaction 
between  strong  mean  shear  and  turbulence  (streaks)  are  the  source  of  thermals  that  penetrate  the  outer 
boundary  layer.  These  coherent  elongated  vertical  motions  interact  with  the  horizontal  mean  wind  to  form 
highly  coherent  secondary  rolls  that  extend  to  the  top  of  the  boundary  layer  and  can  extend  20-40 
boundary  layer  thicknesses  in  the  mean  wind  direction19.  Thus,  there  is  a  direct  coupling  between  the 
near-surface  shear-driven  streaks  and  very  large  eddy  structure  of  the  boundary  layer.  (It  is  possible  that 
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the  mechanisms  that  underlie  the  creation  of  the  longitudinal  convective  rolls  may  be  related  to  the 
mechanisms  underlying  much  weaker  highly  elongated  structures  that  have  been  observed  in  the  log  layer 
of  the  neutral  boundary  layer20.) 

An  important  negative  consequence  of  the  inner-outer  coupling  is  that  the  near-surface  errors  from  the 
overshoot  are  driven  vertically  to  infect  the  entire  boundary  layer.  To  demonstrate  this,  Khanna  & 
Brasseur17  applied  two  SFS  models  to  the  same  LES;  one  produced  a  stronger  overshoot  than  the  other. 
The  LES  with  the  stronger  overshoot  produced  much  stronger  and  coherent  convective  thermals  and 
boundary  layer  rolls  with  much  larger  horizontal  integral  scales  that  persisted  to  the  top  of  the  boundary 
layer.  Furthermore,  these  overly  coherent  thermals  are  spuriously  aligned  with  the  mean  geostrophic 
wind.  This  spuriously  strong  thermal  structure  will  adversely  influence  vertical  transport  to  the  upper 
atmosphere  of  momentum,  thermal  energy,  contaminants,  and  humidity.  Error  in  humidity  predictions 
will  likely  enter  cloud  cover  predictions  and  produce  error  in  solar  radiative  heating  at  the  earth's  surface. 
Incorrect  prediction  of  vertical  transport  of  C02  and  other  greenhouse  gases  may  affect  related  predictions 
of  upper  atmosphere  chemistry. 

The  negative  consequences  of  the  overshoot  in  mean  velocity  gradient  arise  essentially  from  the 
incorrect  prediction  of  Reynolds  stress  anisotropy  near  the  surface.  Juneja  &  Brasseur21  argued  that  the 
incorrect  anisotropy  results  from  a  feedback  interaction  between  the  exaggerated  mean  gradient  and 
Reynolds  stress  production  as  a  consequence  of  inherent  under-resolution  at  the  first  grid  level  that  occurs 
in  LES  of  high  Reynolds  number  turbulent  boundary  layers  when  the  first  grid  level  is  in  the  inertial 
layer.  The  errors  are  exacerbated  by  any  mechanism  that  enhances  vertical  transport,  including 
buoyancy17  and  boundary  layer  separation. 

B.  Previous  studies 

Given  the  fundamental  importance  of  the  overshoot,  there  have  been  a  number  of  studies  that  have 
attempted  either  to  understand  the  cause  of  the  overshoot21,  or  have  attempted  to  remove  it.  Mason  and 
Thomson1,  for  example,  suggested  that  the  overshoot  is  particular  to  eddy  viscosity  closures  where  energy 
is  removed  at  each  point  from  the  resolved  scales  in  contrast  with  the  known  forward/backward  nature  of 
energy  transfer  at  a  point.  To  introduce  “backscatter”  into  the  simulation,  they  added  to  the  resolved 
momentum  equation  a  Langevin-like  stochastic  acceleration  term,  in  addition  to  a  SFS  stress  divergence 
using  the  Smagorinsky  closure  for  SFS  stress.  However,  Mason  &  Thomson1  made  another  significant 
modification — they  reduced  the  eddy  viscosity  near  the  surface  by  making  the  Smagorinsky  length  scale 
proportional  to  z  at  grid  nodes  near  the  surface. 

In  Fig.  1  we  compare  results  from  four  other  groups  of  researchers  between  1994  and  2005  who 
explicitly  addressed  the  problem  of  the  overshoot  in  LES  of  the  ABL.  The  surface  layer  is  shaded  to 
indicate  the  region  over  which  <pm  should  be  predicted  as  a  straight  vertical  line  and  from  which  a  grid- 
independent  prediction  for  mean  velocity  should  emanate.  Fig.  1(a)  shows  predictions  from  Sullivan  et 
aln  who  argued  that,  as  the  surface  is  approached,  the  SFS  stress  closure  should  transition  from  an  eddy 
viscosity  model  appropriate  to  LES  to  a  model  more  appropriate  to  the  Reynolds -averaged  Navier-Stokes 
equation  (RANS).  Although  their  mixed  eddy  viscosity  model  was  purely  dissipative  (no  backscatter),  the 
overshoot  was  significantly  reduced  compared  to  a  pure  SFS  model.  Like  Mason  &  Thomson1,  their 
mixed  model  results  in  a  reduction  in  net  SFS  stress  and  eddy  viscosity  near  the  surface.  Whereas  the 
overshoot  could  be  reduced,  and  perhaps  even  eliminated,  by  the  adjustments  to  the  SFS  model,  a 
dependence  of  (j){z )  on  z  persisted  (Fig  1(a)).  The  Sullivan  et  al.n  model  was  further  modified  by  Ding  et 
al.22.  Recently,  Leveque  et  al 23  modified  the  Smagorinsky  eddy-viscosity  by  subtracting  mean  shear  from 
the  instantaneous  resolved  rate-of-strain  tensor.  Whereas  they  applied  the  model  only  to  low  Reynolds 
number  boundary  layers  with  partially  resolved  viscous  layers,  the  effect  is  again  to  reduce  the  eddy 
viscosity  near  the  surface.  However  the  LES  predictions  failed  to  predict  constant  (j){z)  over  the  entire 
surface  layer. 

Kosovic24  added  additional  nonlinear  terms  in  velocity  gradient  to  the  eddy  viscosity  closure. 
Whereas  major  improvement  in  the  overshoot  in  the  neutral  ABL  was  obtained  (Fig.  1(b)),  grid  resolution 
was  quite  low  and  the  improvement  degraded  at  higher  grid  resolution.  Dynamic  formulations  of  the 
Smagorinsky  eddy  viscosity  closure  have  been  applied  by  Porte-Agel  et  al 25  and  Esau26  to  improve  the 
overshoot  and  capture  the  LOTW.  We  show  in  Fig.  1(c)  results  from  Porte-Agel  et  al 25  who  developed  a 
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scale-dependent  formulation  of  the  dynamic  Smagorinsky  model  that  produces  a  reduction  in  eddy 
viscosity  near  the  surface  and  removes  an  apparent  undershoot  with  the  standard  dynamic  model  (solid 
line).  Whereas  the  modifications  significantly  reduced  the  overshoot,  constant  ^(z)  was  not  obtained 
over  the  surface  layer. 

Chow  et  al .27  combined  a  number  of  modeling  elements,  including  a  dynamic  eddy  viscosity  model28, 
a  “resolvable  subfilter  scale  stress  model”  component29  combined  with  a  deconvolution  procedure30,  and  a 
“canopy  model”31.  As  shown  in  Fig.  1(d),  whereas  over-prediction  of  mean  shear  near  the  surface  could 
be  reduced  with  certain  combinations  of  elements,  it  was  not  clear  which  modeling  elements  were 
responsible,  and  a  robust  grid-independent  solution  was  not  obtained.  A  recent  calculation  by  Drobinski  et 
al ,32,  however,  indicated  that  a  suppression  of  the  overshoot  was  possible  using  a  standard  one-equation 
model  but  with  a  more  refined  grid.  Their  result  will  be  discussed  in  Sects.  VLB  and  VII. B  in  context 
with  the  current  analysis. 

Figure  2  shows  that  what  we  refer  to  here  as  an  “overshoot”  is  described  as  a  “logarithmic  layer 
mismatch”  or  “super  buffer  layer”  in  the  detached-eddy  simulations  of  Nikitin  et  al.14.  Spalart15  and 
Piomelli  and  Balaras16  described  this  as  a  fundamental  unresolved  problem  in  DES.  Whereas  the 
simulations  of  Fig.  1  contain  only  the  inertial  LOTW  layer  due  to  the  presence  of  surface  roughness  and 
the  fluctuating  surface  stress  is  modeled,  the  DES  simulations  used  RANS  to  model  a  viscous  surface 
layer  with  no  slip  at  the  wall  (shown  in  Fig.  2(a))  below  an  inertial-dominated  surface  layer  that  is 
simulated  with  LES,  as  shown  in  Fig.  2(b).  The  logarithmic  layer  mismatch,  shown  by  the  circled  part  of 
the  mean  velocity  profile  Fig.  2(a),  is  shown  in  Fig.  2(b)  to  be  equivalent  to  the  overshoot  phenomenon 
discussed  above. 

Figures  1  and  2  motivate  the  current  research;  we  seek  an  understanding  of  the  essential  mechanisms 
underlying  the  overshoot  that  will  both  eliminate  the  overshoot  and  predict  a  vertical  line  in  ®m  vs.  z 

over  the  entire  surface  layer  with  a  grid-independent  prediction  over  the  entire  boundary  layer. 
Historically,  the  assumption  has  been  that  the  solution  to  the  overshoot  problem  is  strictly  a  closure  issue 
and  all  attempts  to  modify  LES  to  predict  LOTW  have  been  through  adjustments  to  the  SFS  stress  tensor 
model.  Although  there  have  been  significant  advances  made,  the  fundamental  mechanisms  underlying  the 
overshoot  are  not  understood  so  that  a  clear  path  to  robust  grid-independent  LES  that  eliminates  the 
overshoot  and  predicts  LOTW  scaling  has  been  elusive.  We  show  here  that  the  fundamental  issues 
underlying  the  overshoot  and  its  resolution  are  broader  than  the  SFS  model  and  involve  basic 
characteristics  of  the  SFS  stress  closure  integrated  with  the  construction  of  the  LES  grid. 

C.  Useful  clues 

In  the  studies  described  above,  a  number  of  important  observations  have  been  made  that  provide  clues 
to  underlying  issues  and  which  require  explanation: 

1.  The  overshoot  is  influenced  by  the  details  of  the  SFS  model.  This  has  been  discussed  above  in  Sect. 
II. B  (Fig.  1).  Whereas  there  have  been  many  variants  to  the  modeling  process,  all  have  employed  an  eddy 
viscosity  component.  A  common  characteristic  of  the  more  successful  adjustments  to  the  modeling 
process  is  a  reduction  in  eddy  viscosity  near  the  surface  compared  to  unadulterated  models. 

2.  The  overshoot  is  tied  to  the  grid.  Khanna  &  Brasseur18  pointed  out  that  increasing  the  resolution  of  the 
grid,  keeping  all  other  elements  of  the  simulation  unchanged,  does  not  diminish  the  magnitude  of  the 
overshoot,  but  moves  the  peak  in  ®m  closer  to  the  surface  in  proportion  to  the  vertical  grid  spacing  ,  A. 

Similarly,  Spalart15  pointed  out  that  in  DES  “grid  refinement  merely  moves  the  same  amount  of 
[logarithmic  layer]  mismatch  closer  to  the  wall.”  This  dependence  of  the  overshoot  on  grid  resolution  is 
clear,  for  example,  in  Fig.  1(d).  Why  the  location  of  the  overshoot  should  be  proportional  to  Az,  however, 
is  not  understood.  Khanna  &  Brasseur18  pointed  out  that,  because  the  horizontal  integral  scale  of  vertical 
velocity  scales  on  z,  vertical  velocity  is  always  under-resolved  at  the  first  grid  level  independent  of  grid 
resolution  and  that  this  under-resolution  has  negative  consequences  that  must  be  addressed  to  eliminate 
the  overshoot.  Specifically,  they  proposed  that  this  inherent  under-resolution  at  the  first  grid  level 
somehow  ties  the  overshoot  to  the  grid.  The  mechanisms  were  unclear,  but  were  felt  to  be  somehow 
associated  with  the  lack  of  performance  of  SFS  models  when  the  integral  scales  are  under  resolved.  This 
shall  be  discussed  in  Sec.  IV. B. 
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3.  The  prediction  of  mean  velocity  gradient  is  grid-dependent.  A  requirement  for  any  successful 
numerical  simulation  is  grid  independence  in  the  solution  for  mean  variables33.  As  illustrated  in  Fig.  1(d), 
the  solution  for  mean  velocity  gradient  often  does  not  converge  as  the  grid  is  refined;  each  grid  produces  a 
different  solution  not  only  in  the  overshoot  region  but  throughout  the  boundary  layer.  Grid  dependence  in 
the  flow  is  apparent  when  different  solutions  in  the  literature  are  compared  (e.g.,  Andren  et  al .10).  We 
shall  show  in  Sect.  VI.C  that  a  grid-independent  solution  is  only  possible  when  both  the  overshoot  is 
suppressed  and  LOTW  scaling  are  obtained. 

It  is  apparent  from  these  previous  studies  that  although  the  overshoot  is  influenced  by  the  closure  for 
SFS  stress,  the  overshoot  problem  is  only  part  of  the  broader  issue  of  accurately  predicting  the  LOTW, 
and  that  these  are  both  modeling  and  numerical  issues. 


II.  An  Analysis  of  the  fundamental  Nature  of  the  Overshoot: 
the  First  observation 

In  this  section  we  focus  specifically  on  the  overshoot.  A  primary  mechanism  underlying  the  overshoot 
and  its  resolution  can  be  understood  by  comparing  true  inertial-vs. -viscous  scaling  underlying  the 
stationary  fully  developed  smooth-wall  channel  flow  in  the  high  Reynolds  number  limit  with  scaling  of 
large-eddy  simulation  of  the  same  high  Reynolds  number  channel  flow  with  unresolved  viscous  or 
roughness  layer.  We  shall  find  that  we  can  relate  the  true  physics  of  the  channel  flow  to  the  spurious 
physics  of  the  simulated  channel  flow  that  is  model  and  algorithmically  dependent  and  cannot  be  entirely 
eliminated. 


A.  Scaling  high  Reynolds  number  smooth  wall  turbulent  channel  flow 

Consider  a  fully  developed  stationary  smooth  wall  incompressible  channel  flow  at  Reynolds  numbers 
sufficiently  high  to  support  the  classical  LOTW  in  the  surface  layer.  The  local  and  global  mean  axial 
momentum  balances  are,  respectively, 

(2) 

8x  6z  dz  8z  8x  8 


where  T0  =  Ttot  (0)  =  pul  and  8  is  the  half  channel  width,  or  boundary  layer  depth  (the  height  where 
Ttot  =  0 ).  The  second  expression  in  (2)  results  from  a  global  force  balance.  The  coordinates  (v,y,z)  are  in 
the  mean-flow,  spanwise,  and  wall-normal  directions,  respectively.  P(x)  is  the  mean  pressure,  and 
Ttot  =Tt+Tv  is  the  total  mean  stress,  where  Tt  and  Tv  are  Reynolds  stress  and  mean  viscous  stress, 
respectively: 


Tt(z)  =  -p(u'w'),  Tv(z)  =  ju—.  (3) 

oz 

Capital  letters  and  primed  quantities  indicate  mean  flow  and  fluctuating  variables,  respectively,  ( u,v,w ) 
are  the  velocity  components  in  the  (x,y,z)  directions,  (p,ju)  are  density  and  viscosity,  and  the  angle 
brackets  denote  ensemble  averaging. 


Integrating  Eq.  (2)  in  z  and  replacing  Ty  (z)  with 
mean  gradient  is 


uu * 

- ®m(z) ,  the  exact  equation  for  the  normalized 

kz 


<D  m{z)  =  KZ+ 


1-IT(z)-7 

V  o  J 


(4) 


where  Tt+  =Tt  !T0  is  the  ratio  of  the  Reynolds  stress  Tt(z)  to  the  wall  stress  T0,  and  z+  =  zt  £v  is  the 

surface-normal  coordinate  nondimensionalized  by  the  viscous  surface  layer  scale  £v  =  v/u*,  where  v= 

pip.  Equation  (4)  states  that  the  total  stress  Ttot(z ),  split  between  the  viscous  and  inertial  contributions  on 
LHS  and  RHS  of  (4),  decreases  linearly  from  T0  at  the  surface  to  zero  at  the  channel  center-plane.  Near 
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the  surface  (z/ S«  1)  in  the  friction-dominated  layer  (z+  <  5),  the  turbulent  stress  is  negligible  ( 
Tt+  «  1 )  and  Eq.  (4)  is  well  approximated  by 

O m{z)~KZ+  (z  <  5).  (5) 

This  result  suggests  that  <1 >m(z)  exceeds  1  when  z+  >  l/x*~  2.5  (for  0.4),  which  is  well  within  the 
friction-dominated  portion  of  the  surface  layer. 

Equation  (5)  indicates  that  an  overshoot  in  ®m(z)  exists  in  the  smooth-surface  high  Reynolds 

number  channel  flow  in  a  region  bounded  from  below  by  z+  ~  2.5  and  from  above  by  the  lower  margin 
of  the  inertial  LOTW  layer  and  upper  margin  of  the  buffer  layer.  In  Fig.  3  we  show  that  this  is  indeed  the 
case.  We  replot,  in  this  figure,  data  from  direct  numerical  simulations  (DNS)  of  the  smooth-wall  channel 
at  different  Reynolds  numbers  that  has  been  graciously  made  available  to  the  scientific  community  by 
Iwamoto  et  al34,35  at  Rer  =  300,  395,  and  642,  and  at  Rer  =  2003  by  Hoyas  &  Jimenez36,  where 

Rer  =u*S/v .  Figure  3(a)  shows  that  ® m(z+)  exceeds  1  when  z+  exceeds  «2.5  and  peaks  at  z+  ~  10 
with  a  maximum  value  of  about  2.3  independent  of  Reynolds  number.  ®m  approaches  1  at  z+  ~  40-50. 
Comparing  Fig.  3(c)  with  Fig.  3(d),  it  is  significant  that  the  peak  in  Om(z)  occurs  nearly  coincident  with 
the  crossover  between  turbulent  stress  Tt(z)  and  viscous  stress  Tfz).  Both  the  location  of  overshoot  and 
the  location  of  the  crossover  occur  at  «10^v  .  Since  the  position  of  this  overshoot  scales  on  the  viscous 

scale  and  corresponds  to  the  transition  between  the  dominance  of  turbulent  stress  above  and  viscous  stress 
below,  this  overshoot  reflects  the  application  of  an  inertial  length  scale  z  in  the  portion  of  the  surface 
layer  where  the  appropriate  length  scale  is  the  viscous  surface  scale  £v  .  Fig.  3(c)  and  Fig.  3(d)  show  that 
the  overshoot  and  the  crossover  in  Tfz)  and  Tt(z),  move  physically  closer  to  the  surface  with  increasing 
Reynolds  number  without  a  reduction  in  maximum  ®m  . 

B.  LES  of  high  Reynolds  number  turbulent  channel  flow 

We  compare  the  previous  analysis  with  large-eddy  simulation  of  the  same  high  Reynolds  number 
channel  flow  analyzed  in  Sect.  III. A,  but  in  which  either  a  viscous  layer  exists  that  is  fully  unresolved  (£v 

«  Az)  or  the  surface  is  rough  with  fully  unresolved  roughness  scale  ( z0  <^C  A7 ).  The  momentum  equation 
for  the  resolved  velocity  u\  is 

dul  ,  diu'u')  l  drfs 

- 1 - = - ,  (  O) 

dt  dx.  p  dX)  dx. 

where  the  viscous  force  has  been  scaled  out  of  the  momentum  balance  on  account  of  the  high  local 
Reynolds  numbers  on  all  grid  nodes  within  the  computational  domain  (surface  viscous  layers  are 

unresolved  or  nonexistence).  The  subfilter-scale  (SFS)  stress  tensor  t?fs  is  modeled.  We  apply  a 

superscript  r  to  indicate  a  variable  that  is  carried  forward  in  the  simulation  as  a  resolved  variable — that  is, 
after  the  process  of  explicit  and  implicit  filtering  in  the  algorithmic  advancement  of  the  modeled 
discretized  version  of  Eq.  (6).  Explicit  filtering  is  generally  carried  out  algorithmically  as  a  dealiasing 
step1,  typically  in  pseudo-spectral  FES.  Implicit  filtering  arises  from  the  dissipative  nature  of  the  model 

for  t?fs  ,  from  numerical  dissipation  within  the  discretized  version  of  Eq.  (6),  and  from  any  dissipative 
elements  introduced  algorithmically  as  the  discretized  dynamical  system  is  advanced  in  time. 

The  ensemble  mean  of  Eq.  (6)  for  stationary  fully  developed  high  Reynolds  number  channel  flow  is: 


1  Algorithmically,  ( uu )  in  equation  (6)  should  be  written  as  (uuf  to  indicate  the  common  application  of  a 
second  explicit  filter  f  on  the  nonlinear  term. 
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(?) 


dP  dTxiz)  {  8Ts(z)  dT^jz) ,  _itfa  dP  pul 
dx  dz  dz  dz  8x  S’ 

where  T0  =  Ttot  (0)  =  pul .  The  second  equation  is  the  global  force  balance.  Total  stress  Ttot  =TR+TS  is 

now  the  summation  of  a  turbulent  stress  TR  formed  from  the  fluctuating  resolved  velocity  components  and 
the  shear  component  of  the  mean  SFS  stress  tensor,  Ts: 

Tr (z)  =  ~P (“'V' )  .  Ts  (z)  =  ~p{tT  )  •  (8) 

Whereas  the  divergence  of  the  SFS  “stress”  is  physically  an  inertial  contribution  to  the  force  balance,  in 
application  all  SFS  models  are  structured  so  as  to  contain  a  frictional  contribution  to  the  equation  of 
motion.  This  is  explicitly  true  of  eddy  viscosity  models  and  mixed  models  which  include  an  eddy 
viscosity  term. 

In  what  follows,  we  do  not  restrict  our  theoretical  treatment  to  any  particular  SFS  stress  model.  We 
do,  however,  use  eddy  viscosity  closures  for  guidance  and  to  carry  out  numerical  experiments. 

We  seek  a  mechanism  to  extract  the  frictional  component  of  the  complete  modeled  tensor  T*IS  ;  that 
is,  we  seek  an  estimate  of  a  scalar  viscosity  that  extracts  the  part  of  t^fs  that  is  parallel  to  the  resolved 
strain-rate  tensor,  s- ,  in  the  mean.  For  the  purposes  of  determining  the  underlying  causes  of  the 
overshoot  in  the  highly  shear-dominated  part  of  the  boundary  layer,  we  define  a  “LES  viscosity,”  Vies(z)9 
as  the  proportionality  between  —p  (t^s  j  and  2  j  =  dU  /  dz  : 


w 


(closure  independent)  (9) 


where  sr{j  is  the  resolved  strain-rate  tensor.  Furthermore,  since  we  are  focused  on  the  first  few  grid  points 

adjacent  to  the  surface,  we  parameterize  LES  viscosity  with  its  value  at  the  first  grid  level,  VLES  (z, )  ,  and 
we  define  the  normalized  LES  viscosity  as  follows: 


VLES  ~  Vles  (^1 )’ 


les 


0 z)  =  — 


(z) 


LES 


(10) 


We  emphasize  that  the  definitions  in  (9)  and  (10)  do  not  assume  an  eddy  viscosity  model  and  can  be  made 
for  any  closure  of  the  SFS  stress  tensor.  However,  this  estimate  of  frictional  content  does  is  only  valid 
with  strong  mean  shear  at  the  first  grid  level,  and  is  therefore  appropriate  to  the  neutral  boundary  layer, 
and  the  stable  and  moderately  convective  atmospheric  boundary  layers  with  shear-dominated  turbulence 
at  zi. 

dU 

Replacing  Ts  in  Eq.  (7)  with  pyLESvL-(z) - >  integrating  in  z,  and  using  inertial  LOTW  scaling  (Eq. 

dz 

1),  produces  the  following  expression  for  ®m(z)  : 


where 


KZ+  d 

4>„(z)=  “ 


<(z) 


•  A 


1 

V  O 


p  -  Z  P  -  ]jxs  anr|  T+  _  Tr 
zees  —  „  ?  t.  —  ,  ana  i  —  . 


(11) 

(12) 


As  in  the  exact  channel  flow  equations,  =  Ttot(  0)  =  pul ,  so  that  TR(z)  is  the  ratio  of  the  resolved  part 

of  the  Reynolds  shear  stress  to  the  total  wall  shear  stress.  Unlike  the  viscous  layer  where  we  could  argue 
that  the  Reynolds  stress  is  a  small  percentage  of  wall  stress  and  approximate  Eq.  (4)  by  Eq.  (2),  we  cannot 
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a-priori  argue  that  TR  «  T0 .  However,  at  the  first  few  grid  points,  z/ S  « 1  and  v^es{z)  ~  0(1)  ,  so  that 
(11)  can  be  approximated  by 

Om (z)  ~  kz+les  (l-r;)  (first  few  grid  points) .  (13) 

The  magnitude  of  TR  depends  on  the  relative  content  of  resolved  to  SFS  stress  in  total  stress 
Ttot  =Tr+Ts.  When  the  separation  between  resolved  and  subfilter-scales  takes  place  in  the  inertial  range 

of  turbulence  scales,  implying  that  all  integral  scales  are  well  resolved,  then  Ts  «  TR  ,  and  TR  cannot  be 

neglected  in  Eq.  (13).  However,  as  discussed  in  Sect.  II.C,  some  integral  scales  are  inherently  under¬ 
resolved  at  the  first  few  grid  levels  in  high  Reynolds  number  LES  of  wall-bounded  flows.  Thus,  it  may  be 

the  case  that  TR  «  1  close  to  the  surface  so  that  Eq.  (13)  is  approximated  by  ® m(z )  ~  k^les  •  ^  this 
were  the  case,  then  we  would  reach  the  same  conclusion  as  when  Eq.  (4)  was  reduced  to  Eq.  (5)  in  the 
viscous  sublayer  of  the  smooth-wall  channel  flow:  an  overshoot  would  occur  when  zLES  >  1/  ft*~  2.5, 
that  is  when  z  >  2.5A 

~  VLES 

Figure  4  shows  that,  indeed,  Ts  »  TR  near  the  surface  in  a  typical  LES  of  the  neutral  atmospheric 
boundary  layer2  using  the  Smagorinsky  model,  so  that  ®m(z)  «  tczLES  and  an  overshoot  is  produced  at 

ZLEs  >  2.5 .  In  fact,  Fig.  4(a,b)  is  very  similar  to  the  curves  in  Fig.  3(a,b)  of  the  real  overshoot  in  the 

smooth  wall  channel  flow.  The  spurious  overshoot  initiates  between  the  first  and  second  grid  levels  and 
peaks  nearly  coincident  with  the  crossover  between  TR  and  Ts ,  very  similar  to  Fig.  3  where  the  real 
overshoot  is  found  to  be  coincident  with  the  crossover  between  Tt  and  Tv. 

The  overshoot  in  LES  appears  to  arise  from  physics  similar  to  the  true  overshoot  in  smooth  wall 
channel  flow.  However,  in  LES  the  frictional  layer  that  causes  the  overshoot  is  a  numerical  LES  frictional 
layer  near  the  surface  that  arises  from  the  frictional  nature  of  the  modeled  SFS  stress  (and  any  numerical 
and  algorithmic  additions  to  dissipation).  This  conclusion  analyzed  further  in  the  following  sections. 

C.  The  first  criterion 

The  observation  from  Fig.  4  that  the  overshoot  is  associated  with  a  reduction  in  the  resolved  Reynolds 
stress  to  below  the  mean  SFS  stress  suggests  that  the  spurious  frictional  content  of  the  model  for  SFS 
stress  has  introduced  a  spurious  length  scale  £  that  interferes  with  the  inertial  scale  z  that  should 

dominate  the  surface  layer.  The  spurious  interference  of  ■$  with  inertial  scaling  and  the  consequent 

spurious  overshoot  are  essentially  the  same  physics  that  underlie  the  production  of  the  real  overshoot  in 
the  friction-dominated  part  of  the  LOTW  layer  of  the  smooth-wall  channel  flow.  However,  unlike  the  true 
viscous  layer,  which  is  a  necessary  consequence  of  frictional  force  with  no-slip,  the  spurious  LES 
frictional  layer  must  be  controlled  to  eliminate  the  frictionally  induced  spuriously  large  mean  gradients 
near  the  surface.  If  it  were  possible  to  maintain  the  dominance  of  TR  over  Ts  to  the  surface,  then  the 
spurious  frictional  contribution  from  the  SFS  stress  model  would  remain  suppressed.  This  observation 
suggests  a  criterion  for  elimination  of  the  overshoot,  that 

T 

9?  =  —  >  yC  ~  0(1) ,  (closure  independent)  (14) 

ts, 

where  the  subscripts  1  mean  “at  the  first  grid  level.”  The  critical  value  $R*,  to  be  determined 
experimentally,  is  an  order  1  quantity  but  may  depend  on  the  model  for  SFS  stress,  the  lower  boundary 
condition,  the  stability  of  the  ABL,  the  numerical  algorithm,  etc. 


2 


We  describe  the  details  of  our  large-eddy  simulations  of  the  atmospheric  boundary  layer  in  Appendix  B. 
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D.  Understanding  the  ratio  9?  and  the  first  criterion 

On  what  does  9?  depend  and  how  can  it  be  controlled  in  a  large-eddy  simulation?  To  gain  insight  into  the 
nature  of  91  and  find  an  answer  to  this  question,  we  develop  expressions  for  91  based  on  eddy  viscosity 
representations  for  the  deviatoric  part  of  the  SFS  stress,  [t^fs  ]dev  =  rfFS  —\rskFS  /  3]^. : 


[*TLv=-2v^,  v,  =*,*,, 


(15) 


where  £tand  uf  characterize  length  and  energy  at  the  smallest  resolved  scales.  £t  is  often  taken  as  a 
fixed  length  scale  £t  =  Ct A,  where  A  =  (AxAyAz)1/3  is  the  grid  size  ( Ax,  Ay, A7  are  the  grid  spacings  in 
x,y,z)  and  Ct  is  a  model  constant.  Some  models  embed  z-dependence  into  £t  in  the  form  £t(z)  =  Ct£t(z )  , 
where  £t(z)  is  specified  to  decrease  towards  the  surface  and  approaches  A  away  from  the  surface1. 
Dynamic  models  take  £t  =  Ct A,  but  determine  Ct  dynamically  with  the  result  that  Ct  =  Ct(z )  reduces 
near  the  surface25.  Thus,  in  each  case  the  modeled  length  scale  £t(z)  decreases  towards  the  surface.  In 
the  current  work  we  develop  scaling  from  eddy  viscosity  closures  with  constant  £t  =Ct A,  but  shall 
discuss  our  results  in  context  with  z-dependent  £t  in  the  Discussion,  Sect.VII.B. 

The  space-time  variability  in  vt  (x,  t )  is  modeled  through  a  fluctuating  velocity  scale  ut  (x,  t )  ,  given 
for  the  Smagorinsky  closure  by  ut  (x9t)  =  A(2s~s~)1/2 .  In  the  basic  one-equation  model 
ut  x_eq  (x,  t )  =  yje(x9 1 )  ,  where  e(x ,  t )  represents  the  SFS  kinetic  energy,  modeled  through  a  prognostic 
equation37,32.  In  both  models  £t  =  Ct  A .  Here  we  develop  expressions  using  the  Smagorinsky  closure  with 

Ct  replaced  by  C2 ,  as  is  traditional.  However  we  shall  present  results  for  both  the  Smagorinsky  and  one- 
equation  eddy  viscosity  models. 

Independent  of  the  model,  [tx[s  ]dev  =  rffs  and 


where 


Ts,  =-p(rff  ),  =2pvLES{s[,)t  =2P6(t',),(j,'J)|, 

I,  (V.»), 


VLES  =4i  \V,)X  With  4l  = 


(16) 

(17) 


We  have  found  from  LES  of  the  neutral  ABL  that  A  is  typically  ~1 .05.  The  ensemble  mean  of  the  eddy 
viscosity  at  the  first  grid  level  with  the  Smagorinsky  closure  is 


v,)i=C?A2((2^)1/2)i«2C?A2(^)i=C?  A2 


dU  ) 
v  dz 


(18) 


where 


<W'V72(4 3>, 


!  ,  1  («) 
8  (4 


2 

1  J 


is  well  approximated  by  v2  since  the  strain-rate  fluctuations  are  nearly  all  subfilter-scale  (this  is 

especially  true  at  the  first  grid  level  where  the  integral  scales  are  minimally  or  poorly  resolved).  We  write 
the  mean  streamwise  velocity  gradient  in  a  form  appropriate  for  inertial  scaling  in  the  surface  layer3: 


3  When  the  Coriolis  force  is  present,  as  in  LES  of  the  atmospheric  boundary  layer,  the  velocity  at  the  first  grid  level 
is  at  an  angle  6]  to  the  geostrophic  wind  velocity  above  the  boundary  layer.  In  this  case,  equations  (19),  (20),  (21), 
(23),  (24),  (27),  (28)  contain  an  additional  factor  cos#/,  as  given  in  Appendix  A. 
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'dir 
V  dz  J, 


U* 

*1*1 


(19) 


In  Eq.  (19)  kx  is  defined  as  the  value  required  to  make  the  LHS  equal  to  the  RHS  at  the  first  grid  level. 
Only  if  LOTW  is  predicted  by  the  LES,  so  that  k\  is  constant  through  the  surface  layer,  will  k\  be  the 
predicted  value  of  the  von  Karman  constant. 

Inserting  (19)  into  (18)  and  (17)  leads  to  the  following  expression  for  the  LES  viscosity: 

vLES  =  —  D  w* A  ,  where  D  =  C2AR41 3 .  (Smagorinsky)  (20) 

*1 


A„=  A,/A,=A,/A,  is  the  aspect  ratio  of  the  grid  andA7  =  z1  is  the  vertical  grid  spacing.  As  will  be 

discussed  at  length,  Eq.  (20)  indicates  that  the  LES  viscosity  is  proportional  to  the  combination  C^A^4  3, 

and  therefore  is  altered  both  by  the  constant  in  the  SFS  model  and  by  the  aspect  ratio  of  the  grid,  as  well 
as  by  the  vertical  grid  spacing. 

Applying  Kolmogorov  scaling38  and  LOTW  scaling  to  the  one-equation  model  produces  a  result 
similar  to  Eq.  (20),  but  with  Ds  replaced  by  Dk  =  CkAR19  and  with  kx  replaced  by  a  different  order  one 

constant,  /?.  The  main  point  is  that  with  the  eddy  viscosity  closure,  the  LES  viscosity  is  proportional  to  the 
product  of  model  constant  and  grid  aspect  ratio,  each  raised  to  a  power  that  depends  on  the  closure. 

In  order  to  develop  an  expression  for  9?  ,  we  note  that  since  the  total  shear  stress  is  Ttot  =  TR  +  Ts ,  9f 
is  given  by 

9?  _  _m L  _  i  _  ^2 ( Pu *  )  _  I  ?  (closure  independent)  (21) 

Tc  Tn 


where 


V! 


(22) 


In  (21)  and  (22)  Ttot  =  pul  and  Ns  =  S  /  A7  is  the  number  of  grid  points  from  the  surface  to  the  top  of 
the  boundary  layer.  <^2  is  generally  very  close  to  1 . 

Inserting  Eq.  (19)  for  and  Eq.  (20)  for  vLES  into  Ts  =2 pvLES  (s^)  ,  and  then  inserting  this 

result  for  Ts  into  Eq.  (21)  produces  the  following  expression  for  9i  : 

9?  =  - 1 ,  (eddy  viscosity)  (23) 


where  ^  =  is  generally  very  close  to  1.  For  the  Smagorinsky  closure,  /?  =  kx  and 

Dt  Ds  =  C2sAr/?>  ,  whereas  for  the  one-equation  model  Dt  ~^Dk=  CkAR'9  and  [3  is  a  different  order 
one  constant. 

The  eddy  viscosity  closure  was  used  in  the  derivation  of  Eq.  (23)  to  provide  insight  into  the 
mechanisms  underlying  9?  and  how  it  can  be  systematically  adjusted  in  large-eddy  simulation  in  order  to 

move  the  LES  into  the  supercritical  regime  9T  >  9 ?* .  We  learn  that  the  ratio  of  resolved  to  SFS  stress  at 
the  first  grid  level  can  be  increased  either  by  reducing  the  model  constant  or  by  reducing  the  grid  aspect 
ratio  and  that  these  two  changes  act  in  combination  through: 

D,  =  C?ARb ,  (24) 


ll 


where  the  model  constant  Ct  (Eq.  15)  and  AR  enter  in  this  combination  through  the  LES  viscosity,  Eq. 
(20),  and  the  powers  a  and  b  are  model-dependent.  Thus,  LES  viscosity  can  be  decreased  and  91 
increased  by  reducing  either  or  both  C,  and  AR  to  reduce  Dt . 

The  explanation  behind  the  effects  of  model  constant  and  AR  on  91  is  in  the  manner  that  each 
reduction  affects  the  balance  between  resolved  and  SFS  stress  within  the  total  stress  Ttot(z )  that  is  fixed 

by  the  global  momentum  balance.  Reducing  the  model  constant  directly  reduces  the  average  SFS  stress 
Ts.  Reducing  the  aspect  ratio  corresponds  to  an  increase  in  resolution  in  the  horizontal,  which  moves 
vertical  velocity  variance  and  Reynolds  stress  from  subfilter-scales  to  resolved  scales  in  the  horizontal. 
The  consequence  of  both  effects  is  to  reduce  Ts  relative  to  TR,  and  therefore  to  increase  the  ratio 
9 \  =  TrJ  Ts  .  Similarly,  both  effects  reduce  LES  viscosity. 

III.  The  Balance  between  Numerical  Friction  and  Inertia: 
a  Second  Observation 

The  above  suggests  that  the  mechanisms  that  underlie  the  generation  of  a  mean  gradient  overshoot 
and  its  consequences  are  associated  with  numerical  LES  friction  that,  in  the  modeled  dynamical  system, 
forces  a  physical  response  similar  to  that  underlying  the  overshoot  in  Newtonian  turbulent  channel  flow 
(Fig.  3)  that  results  from  molecular  friction.  Correspondingly,  like  the  real  viscous  layer  in  smooth  wall 
channel  flow  that  arises  from  a  change  in  scaling  from  z  to  £v  =  v/u* ,  the  overshoot  in  LES  reflects  a 
transition  in  dominance  from  the  inertial  surface  scale  z  to  a  numerical  LES  viscous  length  scale 
£Vles  =  vLES  /  u*  that  dominates  near  the  surface.  LES  produces  an  overshoot  in  ®m(z)  as  a  result  of  the 

spurious  dominance  of  this  length  scale  £v^  in  a  region  where  the  integral  scales  should  scale  on  z. 

However,  for  the  LOTW  to  exist  in  an  inertia-dominated  surface  layer,  inertial  effects  must  be 
sufficiently  strong  relative  to  viscous  forces,  as  measured  by  the  ratio  of  the  boundary  layer  depth  to  the 
viscous  wall  scale,  the  Reynolds  number  Rer  =  8  /  £ v .  This  suggests  the  existence  of  a  “LES  Reynolds 

number,”  ReL£S  =  5 1  £v^  .  Since  friction  in  the  discretized  LES  dynamical  system  is  at  the  core  of  the 
overshoot,  we  must  consider  all  consequences  of  friction,  including  the  requirement  that  inertial  effects 
dominate  viscous  effects  sufficiently  to  produce  the  LOTW  scaling,  dU  / dz~u*/ z  - 

In  particular,  in  the  smooth  wall  channel  flow  the  inertial  layer  will  only  reflect  the  LOTW  when 
inertia  dominates  the  viscous  force  within  the  surface  layer  sufficiently  that  Rer  exceeds  a  critical  value, 

Re*  .  The  DNS  data  of  Fig.  3  shows  that  the  LOTW  is  not  captured  in  the  viscous  channel  flow  even  at 
Rex  =  2003,  when  the  viscous  layer  (defined  by  the  peak  in  )  is  only  2%  of  the  surface  layer. 

Similarly,  one  can  expect  that  LES  of  the  high  Reynolds  number  boundary  layer  can  only  produce 
LOTW  scaling  when  inertia  in  the  discretized  dynamical  system  dominates  friction  sufficiently  that  the 

LES  Reynolds  number  ReL£5  exceeds  some  critical  value,  Re*EV .  Indeed  one  can  also  postulate  a 
transitional  ReL£5  which  must  be  exceeded  to  support  turbulence  in  the  discretized  LES  dynamical 
system.  In  Fig.  5  we  show  four  large-eddy  simulations  of  the  atmospheric  boundary  layer  at  increasing 
values  of  ReL£5  .  Similar  to  Fig.  3  for  DNS  of  channel  flow  where  the  peak  in  the  overshoot  scales  on  the 

viscous  scale  £v  ,  with  the  exception  of  the  lowest  ReL£5  (thin  solid  curve  with  small  dots),  the  overshoot 
in  LES  scales  on  the  numerical  LES  viscous  scale,  £  (Fig.  5(a)).  In  fact,  both  the  true  channel-flow 
overshoot  and  the  spurious  LES  overshoot  peak  at  10  corresponding  viscous  units.  Thus,  the  LES 
overshoot  moves  closer  to  the  surface  along  with  £  as  ReLES  =  5 1  £  increases  (Fig.  5(c))  and  the 

numerical  LES  viscous  layer  occupies  a  correspondingly  smaller  percentage  of  the  surface  layer. 
Similarly,  the  peak  in  T>m(z)  coincides  with  the  crossover  between  mean  resolved  and  SFS  stress  TR  and 

Ts  (Fig.  5d)  so  that  the  crossover  also  scales  on  £v^  (Fig.  5(b)). 
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The  thin  solid  curve  (with  small  dots)  in  Fig.  5(c)  is  included  to  illustrate  the  consequence  of  a  LES 
Reynolds  number  that  is  too  low  to  support  turbulence.  The  mean  velocity  profile  is  qualitatively  similar 
to  the  parabolic  profile  characteristic  of  laminar  Newtonian  channel  flow.  Thus,  although  the  LES 
equation  contains  no  true  frictional  term,  the  numerical  LES  friction  inherent  in  the  model  and  adjusted 
by  the  grid  as  described  by  vLES  in  Eq.  (20)  can,  like  real  friction,  dampen  inertial  motions  and  prevent 
turbulence  within  the  discretized  dynamical  system  advanced  in  the  LES. 

The  frictional  content  of  the  simulation  embodied  by  vLES  should,  in  principle,  be  extended  to  include 

the  numerical  dissipation  within  the  specific  discretization  that  is  applied  to  advance  the  LES  equations  in 
time  with  a  SFS  model.  The  LES  presented  in  Sect.  VI  apply  the  pseudo-spectral  method  in  the  horizontal 
and  finite  difference  in  the  vertical  on  a  staggered  mesh  and  is  minimally  dissipative.  Significant  frictional 
content  within  the  numerical  algorithm  might  strengthen  the  overshoot  beyond  what  is  described  here. 

A.  The  second  and  third  criteria 

The  first  criterion  discussed  in  Sec.  III.C  is  a  necessary  condition  to  eliminate  the  overshoot  but  is  not 
a  sufficient  condition  to  correctly  predict  LOTW  scaling  in  the  high  Reynolds  number  surface  layer — that 

is,  constant  (j){z )  is  the  surface  layer.  A  second  criterion  is  necessary:  that  ReL£S  exceed  a  critical  Re^ 

to  achieve  LOTW  scaling  in  the  simulated  dynamical  system.  In  the  presence  of  the  overshoot,  one  would 
expect  that  to  predict  the  LOTW  in  the  part  of  the  surface  layer  that  is  not  directly  affected  by  numerical 

LES  friction,  Re^s  would  have  to  achieve  values  in  the  thousands  similar  to  Re*  in  real  frictional 

channel  flow.  However,  unlike  DNS  of  channel  flow  where  the  overshoot  is  real  and  cannot  be 
eliminated,  in  LES  of  high  Reynolds  number  boundary  layers  the  overshoot  and  its  frictional  sources  are 

spurious.  Therefore,  the  critical  Reynolds  number  Re^  required  to  satisfy  the  second  criterion  turns  out 

to  be  much  lower  when  9i  >  9T*  (and  the  overshoot  is  eliminated)  than  when  the  overshoot  is  maintained 
as  in  Fig.  5. 

To  show  this  we  use  Eq.  (19)  to  write: 

Ts,  =  2Pvles  (S13 ),  =  PVLES  •  (closure  independent)  (25) 

'  A  *i*i 

Inserting  Eq.  (25)  into  Eq.  (21)  yields  an  expression  for  ReLES  that  is  valid  independent  of  the  SFS 
model: 

N 

R qles  =  — ~  (9?  + 1) .  (closure  independent)  (26) 

44 

Thus  ReL£5  depends  both  on  9?  and  on  the  vertical  grid  resolution^ .  (Note  that  in  Fig.  5(c)  the  LES 
give  by  the  solid  curve  has  low  ReL£5  because  of  low  vertical  grid  resolution.)  The  critical  LES  Reynolds 
number  Re*£5  is  therefore  associated  not  only  with  the  critical  ratio  of  resolved  to  SFS  stress  %C ,  but 
also  with  a  critical  vertical  resolution,  N*s  where 

*  Nl 

R qLes  ~  “  (9?  + 1)  .  (closure  independent)  (27) 

44 

The  third  criterion,  that  Ns  exceed  a  critical  value  N*s,  follows  from  the  second  criterion, 
ReL^  >  R q*les  when  9T  >  9T* . 

One  way  to  understand  the  requirement  for  a  minimum  vertical  resolution  to  produce  high  accuracy 
LES  of  the  boundary  layer  is  simply  as  a  manifestation  of  the  standard  computational  requirement  that  all 
special  regions  with  their  own  characteristic  dynamics  be  well  resolved  for  accurate  numerical  simulation. 
The  surface  layer  is  an  example  of  a  region  with  special  dynamics  that  requires  good  resolution.  The 
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surface  layer  occupies  -15-20%  of  the  boundary  layer  depth  at  high  Reynolds  numbers.  Resolving  this 
layer  with,  say,  10  grid  points  in  the  vertical  therefore  leads  to  an  estimate  for  N*s  of  -50-65.  Recognizing 

that  9i*  and  Re*EV  are  model-dependent,  we  can  estimate  Re*^  roughly  for  9?*  « 1  to  be 
Re*^  —250—325 .  This  estimate  assumes  that  the  overshoot  has  been  eliminated  and  that  the  LOTW 
has  been  captured  with^  =  /c»0.4.  That  is,  we  assume  that  there  are  no  additional  confounding 
influences  that  cause  the  LES  to  deviate  from  the  LOTW.  (In  Sect.  VII  we  shall  discuss  confounding 
influences.)  This  estimate  for  Re*EV  is  well  over  an  order  of  magnitude  lower  than  what  we  would 

estimate  by  analogy  with  Re*  in  DNS  of  channel  flow  if  the  overshoot  were  retained,  but  confined  to  a 
sufficiently  thin,  numerically  viscous,  layer  adjacent  to  the  surface. 

We  shall  find  (Sect.  VI)  that,  for  our  current  LES  of  the  ABL  with  the  Smagorinsky  eddy  viscosity 
closure,  45-50  is  a  reasonable  estimate  for  N*s  and  350  is  a  reasonable  estimate  for  Re*EV .  The  good 

news  is  that  a  vertical  grid  resolution  «  2N*§  —90-100  is  not  overly  severe  and  doable  on  current 

mainframes.  The  bad  news  is  that  most  calculations  in  the  literature  are  of  LES  with  vertical  resolutions 
below  critical.  Interestingly,  we  shall  show  in  Sect.  V.B  that,  in  addition  to  subcritical  resolution  in  the 
vertical,  there  are  practical  limitations  to  the  maximum  vertical  resolution  in  LES  of  the  high  Reynolds 
number  boundary  layer. 

B.  Understanding  the  LES  Reynolds  number  and  the  three  criteria 

To  develop  greater  insight  into  the  LES  Reynolds  number  and  its  control,  we  evaluate 
ReLES  —  £ yLES  using  the  Smagorinsky  eddy  viscosity  representation  for  the  SFS  stress  as  was  done  to 

understand  9?  in  Sect.  III.D.  There  we  derived  an  expression  for  the  LES  viscosity,  vLES  (Eq.  20).  Using 
this  expression,  the  numerical  LES  viscous  length  scale  is  given  by: 


£ 


VLES 


v_^  =  kDA  , 

w*  J3  ‘  z 


(eddy  viscosity)  (28) 


where  /?  —  kx  for  the  Smagorinsky  closure.  Dividing  8  by  Eq.  (28),  or  replacing  9 \  in  Eq.  (26)  by  Eq. 
(23),  produces  the  following  expression  for  the  LES  Reynolds  number: 

ReLES  =  T  TT  •  (eddy  viscosity)  (29) 

6  D, 

Several  interesting  observations  can  be  extracted  from  Eqs.  (28)  and  (29).  Since  the  overshoot  scales 
on  £  and  peaks  at  10  £  (Fig.  5),  the  observation  made  in  Sect.  II. C  from  previous  studies — that  the 

v LES  v LES 

overshoot  is  tied  to  the  grid — can  now  be  explained.  Equation  (28)  shows  that  if  neither  the  model 
constant  nor  the  grid  aspect  ratio  are  altered  while  the  grid  is  refined,  the  LES  viscous  scale  £  ,  and 

therefore  the  peak  in  ®m  ,  will  move  closer  to  the  surface  in  proportion  to  the  grid  spacing,  Az.  However, 
since  Dt=C?ARb  is  left  unchanged  during  the  grid  refinement,  9i  and  the  magnitude  of  the  overshoot 

will  not  change — the  overshoot  simply  moves  closer  to  the  surface  as  shown  in  Figs.  3  and  5.  In 
particular,  Fig.  5  shows  that  as  the  overshoot  moves  towards  the  surface  in  proportion  with  the  grid 
spacing  at  constant  AR  and  Ch  the  ratio  TR  /  Ts.  =  9T  and  magnitude  of  peak  ®m  remain  unchanged  and 

the  locations  of  peak  ®m  remains  attached  to  the  crossover  between  TR  and  Ts. 

A  second  interpretation  follows  by  replacing  R qles  with  8  /  £*  in  Eq.  (27).  The  criterion 

9t  >  9?*  ~  1  is  then  equivalent  to: 


14 


0.2. 


(closure  independent)  (30) 


VLES  ^  VLES  _  ^2^*1 

Az  Az  ~91*+1 

Equation  (30)  states  that  the  spurious  length  scale  arising  from  friction  within  the  SFS  model  and 

grid  under-resolution  in  a  large-eddy  simulation  must  be  confined  sufficiently  well  within  the  first  grid 
cell  for  numerical  LES  friction  to  not  adversely  affect  the  LES.  Note  that  Eq.  (30)  is  equivalent  to  the 

requirement  that  zf  >  5  :  the  first  grid  level  must  be  at  least  five  times  than  the  spurious  viscous  length 
scale. 

The  inequality  (30)  is  satisfied  when  SH  >  5H* ,  so  that  the  first  and  second  criteria  are  met  when  the 
spurious  viscous  length  scale  is  sufficiently  small  relative  to  the  grid  spacing.  However  (30)  does  not 
guarantee  the  third  criterion,  NS>N*S.  This  is  because  the  overshoot  peaks  at  ^£Vles->  so  that  the 

condition  given  by  (30)  still  allows  partial  resolution  of  the  overshoot  (e.g.,  when  10£ /  A7  ~  2)  .  One 

can  therefore  interpret  the  addition  of  the  third  criterion,  Ns  >  N*s,  as  demanding  that  the  spurious 

frictional  length  scale  £  be  buried  both  sufficiently  far  within  the  first  grid  level  Az  and  within  the 

boundary  layer  8  that  the  grid  is  given  no  opportunity  to  either  create  an  overshoot  or  alter  LOTW  scaling 
through  the  influence  of  friction  within  the  surface  layer. 


IV.  A  Framework  for  High- Accuracy  Large-Eddy  Simulation: 
the  “High- Accuracy  Zone” 

By  comparing  Eq.  (23)  for  3?  with  Eq.  (29)  for  ReL£5  we  learn  that  reducing  the  model  constant 
and/or  the  grid  aspect  ratio  in  the  combination  Dt  =  Cat  Ab  causes  both  91  and  ReL£5  to  increase  (and, 
correspondingly,  £  /  Az  to  decrease).  However,  increasing  only  the  vertical  resolution  increases  the 

LES  Reynolds  number  but  has  no  effect  on  the  ratio  of  resolved  to  SFS  stress,  9?  .  This  observation  leads 
to  the  concept  of  a  “91  —  ReLES  parameter  space”  in  which  high-accuracy  LES  of  the  high  Reynolds 

number  boundary  layer  could  be  developed.  Within  this  framework,  one  can  systematically  adjust  LES  of 
the  boundary  layer  so  that,  in  the  high  Reynolds  number  limit,  the  overshoot  is  suppressed  and  the  LOTW 
is  captured.  The  91  —  ReLES  parameter  space  is  illustrated  in  Fig.  6;  a  large-eddy  simulation  of  the 

boundary  layer  is  identified  as  a  point  on  a  plot  of  9?  against  ReLES  .  In  subsequent  simulations,  the  LES 
is  adjusted  to  move  the  point  within  the  9?  —  ReLES  parameter  space  relative  to  the  critical  parameters 
9?*,  ReLES  and  N*s. 

For  the  LES  to  capture  the  LOTW  while  resolving  the  overshoot,  the  simulation  must  live  in  the 
rectangular  space  5H  >  5H* ,  Re/J?s  >  Re*£S .  We  have  roughly  estimated  91*  ~1  and  Re^5  ~  350. 

However  in  addition  to  the  criteria  9?  >  91*  and  ReL£5  >  Re^ ,  we  have  argued  for  a  third  criterion, 
>  N*s.  9?  and  Re/  f:v  are  linearly  related  by  Eq.  (26): 


9 ?  = 


kNsj 


Re 


LES 


1. 


(closure  independent)  (31) 


Thus,  Ns  enters  in  the  slope  of  91  vs.  ReL£5  .  In  Fig.  6  we  plot  Eq.  (31)  as  a  series  of  lines  with  constant 
slope  ^2kx  /  Ns.  Since  LOTW  is  only  captured  in  the  supercritical  region  of  the  91  —  ReLES  parameter 
space,  in  general  kx  will  vary  from  point  to  point  on  lines  of  constant  slope  <%2kx  / N §.  However  the 
variation  in  kx  is  not  so  great  as  to  obscure  the  strong  inverse  relationship  between  the  slopes  of  the 
91  —  ReLES  lines  and  the  vertical  grid  resolution  Ns-  It  can  be  shown  from  Eq.  (27)  that  the  third  criterion, 
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>  N*s ,  is  met  when  the  simulation  lies  to  the  right  of  the  9i  —  ReLES  line  with  slope  %2kx  /  A/J  that 
passes  through  the  intersection  between  9t  =  9t*  and  ReL£S  =  Re*M  ,  as  illustrated  in  Fig.  6.  We  call  the 
wedge-shaped  region  that  defines  supercritical  LES  satisfying  all  three  criteria  SH  >  SH*  9  R qles  >  Re^  , 

and  Ns  >  N*s  the  “High-Accuracy  Zone,”  or  HAZ.  The  LES  must  reside  within  the  HAZ  to  meet  the 
three  criteria  required  to  both  eliminate  the  overshoot  and  capture  the  LOTW. 

A.  Designing  LES  to  capture  la w-of-th e- wall  scaling 

The  objective  is  to  systematically  move  the  large-eddy  simulation  into  the  HAZ  of  the  9i  —  ReLES 
parameter  space  by  combining  Eq.  (31)  with  knowledge  gained  from  Eqs.  (23)  and  (29).  Although  we 
have  applied  eddy  viscosity  closures  to  gain  insight  into  the  process  of  adjusting  Ns ,  9?  and  ReL£5  to 

systematically  move  the  LES  within  the  9T  —  ReLES  parameter  space,  the  basic  method  can  be  applied 
with  any  SFS  stress  model.  TR  and  7^  (and  therefore  9? ),  and  vLES  (and  therefore  ReLES )  can  be 
quantified  independent  of  SFS  model,  the  only  requirement  being  that  all  contributions  to  the  SFS  stress 
are  included  in  defining  t?fs  before  calculating  7^  .  This  relatively  simple  process  may  be  described  in 
two  basic  steps: 

1.  Adjust,  and  hold  fixed,  the  vertical  resolution  of  the  grid  (Nz)  so  that  when  the  simulation  is  fully 
developed,  Ns  will  exceed  N*s  (typically,  N7  ~  l.5Ns  —2  Ns  to  minimize  the  influence  of  the  upper 

boundary  condition).  We  shall  point  out  in  the  next  section  that  it  is  possible  to  over  resolve  in  the 
vertical. 

2.  Then  systematically  adjust  the  aspect  ratio  of  the  grid  together  with  the  model  constant  to  move  the 
simulation  roughly  along  a  straight  line  in  Fig.  6  from  the  subcritical  region  of  the  9i  —  ReL^ 
parameter  space  into  the  HAZ.  With  eddy  viscosity  closures,  the  model  constant  and  aspect  ratio 
appear  in  the  combination/^  =  C“ARb  .  However  to  move  into  the  HAZ  with  any  other  closure,  one 

would  adjust  the  model  constant  for  that  closure  (to  systematically  reduce  the  SFS  stress)  together 
with  a  systematic  increase  the  horizontal  resolution  of  the  grid  (to  systematically  reduce  the  aspect 
ratio). 

We  presume  that  the  closure  for  SFS  stress  relies  on  a  dissipative  mechanism  to  model  the  net  transfer 
of  resolved  turbulence  energy  to  subfilter  scale  motions.  With  this  methodology  for  designing  LES  one 
can  analyze  systematically  what  works  better  or  worse  depending  on  choice  of  model  type,  model  details, 
model  constant,  grid  resolution,  grid  structure,  algorithm,  geometry,  etc.,  with  some  understanding  of 
underlying  mechanisms.  This  framework  provides  the  LES  community  with  both  physical  understanding 
and  structure  upon  which  a  systematic  procedure  for  LES  design  may  be  based.  Once  the  researcher  has 
become  experienced  with  the  method,  s/he  will  be  able  to  design  high-accuracy  LES  more  rapidly  using 
her/his  favorite  SFS  model,  algorithm,  code,  etc. 

Whereas  the  model  constant  Ct  and  the  model  length  scale  £t  =  Ct A  are  uniform  in  the  above 

discussions,  in  the  general  eddy  viscosity  model  the  length  scale  £t  is  specified  as  varying  with  z.  The 

Mason  &  Thomson1  modification  of  the  length  scale  1 1  near  the  surface  and  the  dynamic  procedure  that 

adjusts  the  Smagorinsky  constant  Cs  with  z  (e.g.,  Porte-Agel  et  al.25)  are  examples.  The  primary  issue  is 
that  the  level  of  eddy  viscosity  be  adjusted  in  concert  with  the  grid  aspect  ratio  (i.e.,  the  horizontal 
resolution  of  vertical  motions)  within  the  first  few  grid  levels  from  the  surface,  where  under-resolution  is 
of  primary  concern  and  mean  SFS  stress  competes  with  resolved  stress. 

B.  Grid-independent  LES  and  practical  limits  on  grid  resolution 

As  discussed  in  Sect.  II.C,  a  problem  with  current  LES  of  the  ABL  is  grid  dependence  in  the  mean 
flow.  We  shall  show  in  the  next  section  that  as  the  simulation  moves  systematically  into  the  HAZ  within 
the  9 1  —  ReL£S  parameter  space,  a  grid-independent  solution  for  the  mean  velocity  is  achieved.  However, 
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one  cannot  move  the  simulation  infinitely  far  into  the  HAZ  along  lines  of  fixed  vertical  grid  resolution 
Ns ,  since  that  would  require  that  either  the  model  constant  be  driven  to  zero  (removing  the  model  from 

the  dynamical  system)  or  the  grid  aspect  ratio  would  be  taken  to  zero  (creating  infinitesimally  thin  grid 
cells  and  infinite  computational  expense).  Either  of  these  limits  will  cause  numerical  problems  and 
simulation  error  regardless  of  computational  expense.  Thus,  for  both  accuracy  and  practical  reasons,  the 
optimal  location  for  the  simulation  within  the  HAZ  is  near  the  apex  of  the  wedge  in  Fig.  6  where 

9?  >  9?  ,  Rg/  7,  S’  >  R ^les  ’  an(^  >  Ns . 

Similarly  there  are  practical  limits  on  vertical  resolution  that,  surprisingly,  confine  the  growth  of  LES 
grids  on  any  given  computer.  Although  a  minimum  vertical  resolution  is  required  to  move  the  simulation 
into  the  HAZ,  progressive  increases  in  vertical  resolution  will  force  the  simulation  onto  lines  in  the 
9T  —  ReLES  parameter  space  that  have  progressively  smaller  slopes,  as  illustrated  by  the  dashed  line  in 
Fig.  6.  Increases  in  vertical  grid  resolution  in  the  absence  of  other  adjustments  will  increase  grid  aspect 
ratio,  driving  91  to  sub  critical  values  and  restoring  the  overshoot  along  with  its  related  errors.  With  eddy 
viscosity  models,  for  example,  in  order  that  91-91*  as  Ns  is  increased,  Dt  =  C“AbR  must  be  held 

constant  (Eq.  23).  However,  since  the  minimum  model  constant  and  maximum  grid  aspect  ratio  are 
bounded,  grid  aspect  ratio  must,  at  some  point,  be  maintained  roughly  fixed  with  increasing  vertical 
resolution,  and  horizontal  grid  resolution  must  increase  proportionally  with  Ns .  The  number  of  grid 

points  will  therefore  increase  approximately  as  N3§ ,  severely  limiting  the  vertical  resolution  to  modest 

values.  This  dilemma  is  reminiscent  of  direct  numerical  simulation  where  the  highest  Reynolds  number 
that  can  be  simulated  accurately  grows  slowly  with  increasing  computer  size  due  to  the  rapid  increase  in 
resolution  requirements  with  increasing  Reynolds  number. 

V.  Numerical  Experiments 

To  evaluate  the  theory  and  further  explore  the  application  of  the  9T  —  ReLES  framework  to  the 

development  of  wall  bounded  LES,  we  have  carried  out  over  110  large-eddy  simulations  of  the  neutral 
shear-driven  atmospheric  boundary  layer  capped  with  an  inversion  layer  to  suppress  boundary  layer 
growth  and  produce  a  quasi-stationary  long-time  solution.  The  Coriolis  force  is  included  at  a  relatively 
high  level  to  reduce  the  time  to  reach  quasi  stationary,  so  the  mean  wind  is  skewed  relative  to  the 
geostrophic  wind  (x  direction)  at  the  first  grid  level  (see  Appendix  A).  The  code  is  pseudo-spectral  in  the 
horizontal  and  finite  difference  in  the  vertical,  so  numerical  dissipation  is  minimal.  In  the  horizontal 
statistically  homogeneous  directions  we  apply  periodic  boundary  conditions;  in  the  vertical  we  apply  the 
boundary  conditions  as  described  in  Moeng37  and  Sullivan  et  al39.  We  report  here  on  simulations  with  the 
Smagorinsky  closure  and  uniform  grid  spacing.  In  particular,  we  apply  the  nonlinear  Moeng37  model  for 
total  fluctuating  shear  stress  at  the  lower  surface  and  the  friction  velocity  is  made  proportional  to  the 
mean  wind  at  the  first  grid  level  with  a  proportionality  constant  that  can  be  related  to  the  surface 
roughness  length  scale,  z0 .  A  summary  of  the  numerical  algorithm  and  simulations  is  given  in  Appendix 

37  11  39  18 

B  and  in  several  publications  ’  ’  ’  . 

It  should  be  noted  that,  in  our  code,  dealiasing  in  the  horizontal  directions  is  carried  out  by  padding 
rather  than  truncation.  These  are  equivalent  except  in  the  interpretation  of  the  grid  and  grid  length  scale: 
the  grid  resolutions  and  aspect  ratios  quoted  in  the  figures  are  before  padding  and  the  grid  length  scale 

A  =  (AxAyA7)l/3  used  in  the  in  the  eddy  viscosity  model  (Eq.  15)  is  also  based  on  the  pre -padded  grid.  In 

all  plots  we  define  the  boundary  layer  thickness  5  at  the  height  where  the  mean  velocity  gradient  crosses 
zero.  Whereas  the  number  of  grid  points  in  the  vertical,  Nz,  is  defined  before  the  simulation,  the  number 
of  grid  points  within  the  boundary  layer,  N&  is  determined  after  a  simulation  is  analyzed  in  the  quasi- 
stationary  state.  We  were  careful  to  determine  the  quasi-stationary  state  consistently  in  all  simulations 
(see  Appendix  B).  In  all  plots  of  ®m(z)  =  we  use  a  value  of  0.40  for  k. 
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A.  The  9T  —  ReLES  parameter  space 

In  Fig.  7  we  identify  in  the  9f  —  ReLES  parameter  space  all  the  simulations  carried  out  for  the  current 
study.  As  will  be  discussed  below,  from  the  predictions  of  ®m(z)  we  estimated  the  critical  lines 
91  =  9?*,  ReLES  =  Re*£S, ,  and  Ns  =  N*s  drawn  on  the  figures.  In  Fig.  7(a)  different  symbols  are  used 
according  to  the  number  of  grid  points  in  the  vertical,  Nz  («  1.5  Ns).  These  show  approximate 
correspondence  to  the  straight  line  representations  of  Eq.  (31)  in  Fig.  6  at  constant  slope,  {£,2^\  /  Ns)  . 
Note  that  the  points  in  Fig.  7(a)  become  noticeably  more  linear  as  the  simulations  enter  the  triangular 
HAZ  region  of  the  9f  —  ReLES  parameter  space  and  kx  better  represents  the  von  Karman  constant.  To 

demonstrate,  as  per  Eq.  (23),  that  9?  increases  with  decreasing  Ds  =  C2AR413,  in  Fig.  7(b)  different 
symbols  are  used  according  to  C2AR4/3  creating  horizontal  bands  of  distinct  symbol  type  consistent  with 

increasing  9?.  Note  that  the  bands  also  shift  to  higher  ReL£5  with  increasing  9?  and  C2AR4/3  as  per 
Eqs.  (26)  and  (29). 

B.  The  “High  Accuracy  Zone”  (HAZ) 

To  show  the  transition  from  subcritical  regions  of  the  9?  —  ReLES  parameter  space  to  the  supercritical 
HAZ,  consider  the  16  simulations  on  the  four  paths  shown  in  Fig.  8  with  the  point  symbols.  Each  path  has 
fixed  vertical  resolution,  Nz,  and  progresses  from  lower  to  higher  9?  and  ReL£5  along  lines  of  constant 

Nz ,  roughly  lines  of  constant  slope  (£,2kx  /  Ns)  in  Eq.  (31). 

In  Fig.  9  we  plot  ®m(z)  against  zJ8  for  each  of  the  groups  of  4  simulations  at  fixed  Nz  in  Fig.  8.  In 

Sect.  IV.A  we  argued  for  the  existence  of  a  critical  vertical  resolution  N*s  below  which  the  surface  layer 

is  insufficiently  resolved  and  the  LOTW  can  therefore  not  be  properly  captured.  The  four  groups  of 
curves  in  Fig.  9  show  the  consequences  of  poor  resolution  of  the  surface  layer  in  the  vertical  and  provides 

an  estimate  for  N*s .  In  Fig.  9(a)  the  surface  layer  is  resolved  with,  at  best,  3  grid  points.  This  clearly 

insufficient  vertical  resolution  is  associated  with  ®m(z)  curves  that  incorrectly  bend  to  the  left  from  the 

surface  with  increasing  z\  this  occurs  even  when  9i  exceeds  1 .  Furthermore,  the  under-resolution  of  the 
surface  layer  in  the  vertical  affects  the  mean  velocity  gradient  and  mean  velocity  profiles  throughout  the 
boundary  layer;  the  mean  gradient  decreases  too  rapidly  in  z-  It  is  not  until  Fig.  9(c),  Nz  =  96,  that  this 
spurious  drop  in  ®m  with  z  transitions  to  the  profile  ®m(z)  extending  vertically  from  the  surface  as  is 
required  by  LOTW.  The  number  of  grid  points  in  the  surface  layer  (up  to  zJ8=  0.2)  is  8-  9,  suggesting 
that  this  is  the  critical  resolution  for  the  surface  layer  and  that  N*s  for  these  large-eddy  simulations  is  ~ 

40-45.  We  have  drawn  a  line  at  constant  slope  k!  N*s  =  0.4/ 45  in  Figs.  7  and  8  showing  how  paths  (a) 
and  (b)  remain  outside  the  HAZ  and  paths  (c)  and  (d)  enter  the  HAZ  at  higher  9?  and  ReL£5  . 

By  comparing  the  groups  of  ®m  curves  (c)  and  (d)  in  Fig.  9,  we  note  that  as  the  simulations  progress 
to  higher  91  and  ReL£5  along  lines  of  constant  Nz,  in  both  groups  the  final  two  curves  approach  grid 
independence  and  the  LOTW  is  captured  as  the  overshoot  is  suppressed.  To  study  the  transition  in  the 
surface  layer  in  detail,  and  to  estimate  9f  and  Re*EV ,  we  compare  the  six  simulations  shown  in  Fig.  8 
with  the  open  circles  which  lie  along  each  of  the  two  paths,  Nz  =  96  and  Nz  =  128.  The  variations  in 
®m(z)  in  the  lower  40%  of  the  boundary  layer  for  each  point  along  the  two  paths  are  shown  separately 
in  Fig.  10.  Comparing  the  upper  and  lower  sets  of  figures,  we  observe  nearly  identical  transitions  along 
each  of  the  two  paths  into  the  HAZ.  At  the  lowest  91  and  ReL£5  an  overshoot  obscures  the  surface  layer 

and  LOTW  is  not  predicted  by  the  LES.  As  9 1  and  ReL£5  increase  (by  decreasing  C2A413 ),  the  same 
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transition  in  0(z)  takes  place  along  each  path:  the  overshoot  progressively  decreases  and  a  region  of 
constant  O m(z)  progressively  strengthens  until  a  well-defined  region  of  constant  ®m(z)  appears  in  the 

surface  layer  that  is  maintained  when  9i*  and  Re^  exceed  about  0.85  and  350,  respectively,  using  the 
Smagorinsky  closure.  (The  oscillations  near  the  surface  will  be  discussed  in  the  next  section.)  The 
estimates  for  the  lines  that  define  the  HAZ  from  the  critical  Smagorinsky  model  parameters  «  0.85 , 
R e*LES  ~  350  ,  and  N*s  ~  45  are  drawn  in  figures  7,  8,  and  1 1(a). 

Figure  10  shows  that  as  Nz,  Cs  and  AR  are  adjusted  so  as  to  progressively  move  the  large-eddy 
simulation  from  the  subcritical  part  of  the  parameter  space  into  the  High-Accuracy  Zone,  the  mean 
velocity  gradient  approaches  a  fixed  point  independent  of  the  path  taken  into  the  HAZ.  In  Fig.  11  we  plot 
normalized  mean  velocity  gradient  distributions  for  five  of  our  computed  large-eddy  simulations  within 

the  HAZ  at  different  vertical  grid  resolutions  Ns  and  combinations,  C2AR413  .  Fig.  11(a)  shows  the 
locations  of  these  five  simulations  on  the  9T  —  ReLES  parameter  space  (all  but  the  +  symbol).  Fig.  11(b) 

shows  that  these  simulations  within  the  HAZ  are  approximately  grid-independent  over  the  boundary  layer 
depth.  In  Fig.  1 1(c)  we  expand  the  surface  layer  to  show  how  the  overshoot  is  suppressed  and  the  LOTW 
is  captured  by  the  simulation  (all  but  the  dashed  line). 

In  Fig.  7(a)  we  also  show  simulations  from  the  literature.  All  but  the  simulation  by  Drobinski  et  al32 
are  well  outside  the  HAZ.  In  two  simulations,  one  from  Andren  et  al.10  and  one  from  Porte-Agel  et  al.25, 
was  sufficiently  high  to  remove  the  overshoot  but  the  vertical  resolution  of  the  grid  was  insufficient 
for  the  LES  to  produce  LOTW  scaling. 

In  Fig.  1 1(a)  the  Drobinski  et  al 32  simulation  is  indicated  with  the  plus  symbol  and  in  Fig.  1 1(c)  with 
the  dashed  line.  Drobinski  et  al 32  used  a  different  numerical  algorithm,  a  one-equation  eddy  viscosity 

model.  Their  model  for  surface  shear  stress  was  not  indicated.  Whereas  their  critical  parameters  9J*  and 
R q*les  are  not  necessarily  the  same  (since  they  used  a  different  SFS  stress  closure),  their  simulation 
appears  to  be  well  within  the  HAZ  (Fig.  11a),  their  overshoot  is  suppressed  and  they  have  captured 
LOTW  scaling  (Fig.  1  lc).  The  Drobinski  et  al 32  simulation  predicts  a  different  numerical  value  for  (f>(z) 
in  the  surface  layer  than  do  our  simulations  with  the  classical  Smagorinsky  closure  and  Moeng37  surface 
stress  model.  Whereas  the  Smagorinsky  model  simulations  in  Fig.  1 1  predicted  a  von  Karman  constant  k 
«  0.33,  the  Drobinski  et  al32  simulation  produces  a  0.35.  Factors  that  affect  the  prediction  of  the  von 
Karman  constant  are  discussed  in  section  VII. C. 

It  turns  out  that  the  von  Karman  constant  issue  is  associated  with  the  development  of  oscillations  in 
®m(z)  at  grid  cells  adjacent  to  the  surface  (Fig.  10).  The  oscillations  grow  as  the  simulations  move 

farther  into  the  HAZ  along  a  line  of  constant  Nz.  As  will  be  discussed  in  Sect.  VII,  the  oscillations 
originate  in  the  spurious  nature  of  the  lower  boundary  condition  for  fluctuating  surface  shear  stress.  This 
observation  is  interesting,  in  part,  because  it  has  been  previously  reported  that  LES  of  the  ABL  is 
insensitive  to  the  details  of  the  lower  stress  boundary  condition11,40.  These  simulations,  however,  were 
outside  the  HAZ  and  in  the  presence  of  the  overshoot  and  its  frictional  source. 
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FIG.  1.  Examples  of  the  overshoot  in  mean  shear  from  previous  LES  studies:  (a)  Sullivan  et  al.n. 
(b)  Kosovic24.  (c)  Porte- Agel  et  al25,  and  (d)  Chow  et  al.21 .  In  (a)  z,  is  the  ABL  depth  defined  in 
the  traditional  manner  as  the  height  where  vertical  turbulent  heat  flux  is  a  minimum.  In  (b)  z  is 
scaled  on  ujf ,  where /is  the  Coriolis  parameter  (the  angular  velocity  at  the  earth's  surface). 
The  simulations  in  (c)  and  (d)  are  pseudo- ABL/channel-flow  simulations  which  do  not  contain  a 
capping  inversion  and  instead  apply  fixed  horizontal  velocity  at  z  =  H.  k  =  0.4  was  assumed  in 
forming  Om . 
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FIG.  2.  The  “logarithmic  layer  mismatch”  in  the  smooth-wall  turbulent  channel  flow  detached- 
eddy  simulations  of  Nikitin  et  a/.14,  discussed  by  Spalart15.  Curve  A4  is  the  simulation  in  (b); 
Rer  =20,000.  (a)  U+  =U / u  plotted  against  z+  =  zl  (,v.  (b)  The  “log  layer  mismatch”  shown 
by  the  upper  oval  in  (a)  is  shown  here  to  be  an  overshoot  in  normalized  mean  shear  Om ,  negating 
LOTW  scaling  over  the  surface  layer.  Figure  (b)  was  kindly  generated  by  Spalart  (private 
communication),  k  is  assumed  to  be  0.4. 
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FIG.  3.  DNS  of  the  smooth-wall  channel  flow  showing  the  physically  real  overshoot  in  the 
viscous  region,  (a)  Normalized  mean  shear  <J>m  vs.  z+ ■  (b)  Wall-normalized  Reynolds  shear 
stress  Tt+  (filled  symbols)  and  mean  viscous  shear  stress,  T.'  (open  symbols)  vs.  z+ .  The  dotted, 
dashed  and  thin  straight  lines  are  the  total  shear  stress  profiles  for  Reynolds  number  of  300,  395, 
and  642,  respectively,  (c)  Omvs.  z/S,  where  5  is  the  half-channel  height,  (d)  T+  (filled 
symbols)  and  T)  (open  symbols)  vs.  z/S .  Total  shear  stress  profiles  for  three  different  Reynolds 
number  collapse  on  to  a  single  linear  line,  k  =  0.4  is  assumed  in  forming  Om  .  Rer  =  300,  395 
and  642  DNS  data  are  from  Iwamoto  et  al.34  and  Iwamoto35;  In  (c)  we  have  added  the  Rer  = 
2003  DNS  data  from  Hoyas  &  Jimenez  .  The  horizontal  dashed  lines  in  (c)  and  (d)  indicated  the 
upper  margin  of  the  surface  layer  where  LOTW  should  be  predicted. 
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FIG.  4.  Overshoot  in  LES  of  the  neutral  ABL  using  the  Smagorinsky  model  ( Cs  =  0.2 )  and  a 
128x128x128  grid,  (a)  ®m ,  vs.  z+LES .  (b)  Wall  normalized  mean  resolved  and  SFS  shear  stress, 
rs+ ,  and  Ts+  plotted  against  z+LES ,  and  its  sum.  k  =  0.4  is  assumed  in  forming  Om .  k  =  0.4  is 
assumed  in  forming  <E>m  . 
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FIG.  5.  Overshoots  in  LES  of  the  neutral  ABL.  (a)  ®m  vs.  z^es  •  03)  Normalized  resolved 
Reynolds  stress  TR  (filled  symbols)  and  mean  SFS  shear  stress  7)+  (open  symbols)  vs.  zEES  . 
Dotted,  dashed  and  thin  lines  are  the  sum  of  resolved  and  SFS  stress,  from  low  to  high  LES 
Reynolds  number,  (c)  <t>m  vs.  z  /  S ,  where  8  is  defined  as  the  height  where  ®m  =  0 .  (d)  TR 
(filled  symbols)  and  TR  (open  symbols)  vs.  z  /  8 .  The  LES  Reynolds  numbers  of  the  simulations 
are  shown  in  (a).  In  order  of  ReL£S ,  the  Smagorinsky  constants  and  grids  were:  ( 
Cs=  0.1,42x42x96),  (^=0.2,192x192x128),  and  (C,  =0.1, 128x128x256).  The  thin 
black  line  in  (c)  is  a  simulation  with  such  low  Re/fS  that  turbulence  is  barely  sustained  ( 
Cs  =  0.2,  42x42x32  and  ReL£V  =  38.).  The  horizontal  dashed  lines  in  (c)  and  (d)  indicated  the 
upper  margin  of  the  surface  layer  where  LOTW  should  be  predicted,  k  =  0.4  is  assumed  in 
forming  ®m . 
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FIG.  7.  The  simulations  carried  out  for  this  study  shown  on  the  '.R  -  ReLES  parameter  space,  (a) 
The  LES  grouped  according  to  Nz,  where  N_  =  32(«),  64("),  96(A),  128(-fc),  160(^)  and 
256(+).  Also  included  are  LES  from  several  previous  studies  by  Andren  et  al.  1994  (O), 
Sullivan  et  al.  1 1  (A),  Porte-Agel  et  al.25  (□),  Chow  et  al.  2005  27  (0)  and  Drobinski  et  al.32  (V). 
(b)  The  LES  grouped  according  to  D  =  C2AR4/3 . 
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FIG.  8.  The  systematic  variations  on  the  9!  -  Re/ES.  parameter  space  for  the  simulations  plotted  in 
Figs.  9  and  10.  Fig.  9(a)  shows  four  simulations  with  the  same  vertical  resolution,  Nz  =  32,  but 
decreasing  Cs2AR4'3  represented  by  +,  x,  ■.  A.  The  vertical  resolutions  in  Figs.  9(b,c,d)  are 
N_  =64,  96,  128,  respectively.  The  open  circles  show  the  locations  of  the  6  simulations  in  Figs. 
9(a,b).  The  vertical  resolutions  in  Figs.  10(a,b)  are  N-  =  96,  128,  respectively. 
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FIG.  9.  Effect  of  vertical  grid  resolution  on  the  LES  predictions  of  mean  shear.  Refer  to  Fig.  8 
for  the  locations  of  the  individual  simulations  on  the  S.H  -  ReL£S  parameter  space.  9?  and  ReLES 
progressively  increase  in  each  figure  in  this  order:  +,  x,  ■,  A.  (a)  N_  =  32,  Nx  =Ny  increased 
from  32(+)  to  128(A).  (b)  Nz  =64,  Nx  =Ny  increased  from  64(+)  to  240(A).  (c)  Nz  =96, 
Nx  =  Ny  increased  from  64(+)  to  288(A).  (d)  N_  =128,  Nx=Ny  increased  from  64(+)  to 
360(A).  k  =  0.4  is  assumed  in  forming  Om .  S  is  defined  as  the  height  where  <E>m  =  0 .  The 
dashed  horizontal  lines  indicated  the  upper  margin  of  the  surface  layer. 
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FIG.  10.  Change  of  structure  of  ®m(z)  as  the  LES  move  towards  and  into  the  HAZ.  Refer  to  Fig. 
8  for  the  locations  of  the  individual  simulations  on  the  91  -  Rei£S,  parameter  space;  91  and  ReLES 
progressively  increase  in  each  panel  from  left  to  right.  Top  panel:  LES  with  Nz  =96.  Bottom 
panel:  LES  with  Nz  =  128 .  ReLES  and  91  are  given  on  top  of  each  figure,  as  ( ReL£S ,  91 ).  k  =  0.4  is 
assumed  in  forming  <E>m  .  5  is  defined  as  the  height  where  ®m  =  0 .  The  horizontal  dashed  lines 
indicate  roughly  the  upper-most  margin  of  the  surface  layer. 
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FIG.  11.  Five  simulations  in  the  HAZ  at  different  Ns  and  Ds  =  C2SA4'3 ,  demonstrating  the 
convergence  of  the  predictions  of  d>m(z)  to  a  relatively  grid-independent  solution  without  the 
overshoot  and  capturing  the  LOTW.  Figure  (a)  gives  the  locations  of  the  5  simulations  in  HAZ. 
Om(z)  from  Drobinski  et  al.32  are  given  in  (c)  with  the  dashed  curve  and  in  (a)  with  the  plus 
symbol  (+).  k  =  0.4  is  assumed  in  forming  Om.  c>  is  defined  as  the  height  where  <E>m  =  0  . 
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